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ABSTRACT 


Three dimensional Inelastic response of laminated composites to 
thermal and mechanical loading Is formulated and analyzed using the 
finite element method. Individual laminae are modeled as homogeneous 
and transversely Isotropic, with material nonllnearlties Introduced via 
a Hill -type yield criterion and an Incremental plasticity approach. 
Nonlinear isotropic hardening is described by Ramberg-Osgood representa- 
tions of uniaxial stress-strain data. A nonlinear hardening coefficient 
is derived based on extrapolation of plastic straining during the 
incremental loading process. Constitutive equations for thermal effects 
include either constant or temperature dependent expansion coefficients. 

The formulations are used as the basis for development of a com- 
puter program designated NALCOM (Nonlinear Analysis of Laminated 
COMposites). The program uses a fully three-dimensional isoparametric 
finite element with 24 nodes and 72 degrees of freedom. An incremental 
solution is performed with nonlinearities introduced as pseudoloads 
computed from initial strains. Equilibrium iteration may be performed 
at every step. Elastic and elastic-plastic response of boron/epoxy and 
graphite/epoxy and problems of curing [0/90] s Gr/Ep laminates with and 
without circular holes are analyzed. Mechanical loading of [±45] s Gr/Ep 
laminates is modeled and symmetry conditions which exist in angle-ply 
laminates are discussed. Results are compared to experiments and other 
analytical models when possible. All models are seen to agree reason- 
ably well with experimental results for off-axis tensile coupons. The 
laminate analyses show the three-dimensional effects which are present 
near holes and free corners. Symmetry conditions often used for 
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angle-ply laminates are shown to be substantially In error for thick 
laminates* becoming better approximations as the thickness decreases. 

An appendix discusses the computer Implementation of the model and 
Includes storage allocations, a flow chart, and a user's guide. 
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Chapter 1 
INTRODUCTION 

Since the Initial work of Smith [1] for plywood, considerable 
research has been conducted In the area of stress analysis of laminated 
orthotropic materials. With the Introduction of fiberglass and the 
advanced (e.g. boron and graphite) and Intermediate (e.g. Kevlar) 
composites, efforts were accelerated due to their potential In aero- 
space applications. Composite materials exhibit higher stiffness to 
weight ratios than their Isotropic counterparts. They can be mad? 
highly corrosion resistant. Compos Ite materials are not, however, 
free of problems. Due to the stiffness mismatch of the orthotropic 
layers bonded together with different fiber orientations. Internal 
stresses, known as Interlaminar stresses, may be developed In the 
vicinity of free edges. For certain stacking sequences and loading 
conditions Interlaminar effects may be the dominant factor In 
structural design. 

Analyses of varying degrees of complexity have been applied to 
laminated composites. The simplest techniques have limited areas of 
applicability, but perform well when properly applied. The more 
complex techniques are computer based, often allowing for detailed 
material models, geometry, and loading conditions. They also *-e 
accurate when properly applied. The increased complexity is expensive, 
both In terms of machine time and man time for setup and later for 
interpretation of results. 

Considerable progress has been made since Smith's efforts, but 
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several areas In the analysis of composites seem to be lacking. 
Composites are known to exhibit significant nonlinearity In stress- 
strain behavior, even at relatively low strains. The nonlinearity 
Is not Isotropic, but varies with direction, as do the elastic 
properties. Models for such elastic-plastic behavior of orthotropic 
materials are not well developed. Two dimensional models for 
nonlinear material behavior have been developed, as have three 
dimensional linear elastic models. Time dependent straining (creep) 
phenomena have been modeled to a limited extent, but more research 
is needed. 

The present study presents a model for the three dimensional 
elastic-plastic stress analysis of layered orthotropic materials. 
Varying degrees of thermal nonlinearity are considered. A computer 
program developed for implementation of the model is validated on 
problems involving nonlinear stress analysis of laminated composites 
including finite width laminates with and without central holes 
subjected to a variety of thermal and mechanical loads. Particular 
attention is given to the interlaminar stresses near free edges. 


Chapter 2 

COMPOSITE MATERIALS LITERATURE REVIEW 

The early work In the analysis of composite laminates was prompted 
by the need to design optimal plywood structures. In 1953, Smith [1] 
presented a means of computing the effective shear modulus for plywood. 
A state of plane stress was assumed. Relssner and Stavsky [2] 
presented a means of Including coupling between In-plane stretching 
and transverse bending. This work presents the basis of the so-called 
"lamination theory" for analysis of materials composed of laminated 
orthotropic plies. Ashton and Whitney [3] present a detailed 
description of lamination theory and Its extension to dynamic and 
stability regimes. Jones [4] also discusses lamination theory in 
detail . 

Considerable effort has gone Into prediction of failure of 
laminated composites. Failure models range from very simple maximum 
stress or strain models to sophisticated models which allow progressive 
failure of individual plies and concurrent assumption of load by 
unfailed plies. Failure prediction is beyond the scope of this study, 
and will not be discussed further. For more information on failure 
criteria of laminates, the reader is directed to the book edited by 
Herakovlch [5]. 

The basic areas of composites discussed herein are divided between 
lamination theory studies and a more detailed group of solutions 
referred to here as full-field solutions. The nature of the two 
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areas and rationale for the nomenclature will be evident as the areas 
are discussed. 

As previously mentioned, the Initial work on lamination theory was 
presented by Relssner and Stavsky [2]. The basic assumptions of 
classical lamination theory (CLT) are that the laminate Is In a state 
of plane stress, and the Individual laminae act Independently. Thus, 
Integrals through the thickness of a laminate are replaced by 
Integrals over the Individual laminae with a summation over the 
laminate. Two significant effects arise from this procedure. First, 
all interlaminar effects are omitted. That Is, all forces which exist 
between the laminae are ignored. Such forces may be shown by more 
sophisticated numerical analyses to be significant. Second, the 
effects of stacking sequence on stress distributions are lost. These 
effects have also been shown to be significant near free edges. 

In defense of CLT, it should be noted that away from free edges 
the CLT solutions are recovered exactly even by the most sophisticated 
techniques. The area of non-applicability of CLT Is confined to a 
region which extends inward from the edge to a distance approximately 
equal to the laminate thickness. This region Is normally referred to 
as the boundary layer. Pagano [6,7] has used a modified lamination 
theory to model the boundary layer. Stresses are shown to change 
rapidly with spatial position and grow very large near free edges. 
Several investigators [8,9] have postulated the existence of stress 
singularities In the boundary layer for linear elastic material 
behavior. In reality, singular stresses do not exist, but are 
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relieved by plasticity or failure. Thus, accurate modeling of free edge 
effects requires an Inelastic capability, and possibly the Inclusion of 
failure mechanisms. 

The early lamination theory work was presented for linear elastic 
materials with no thermal effects. Tsai [10] presented a means of 
Including thermal loading based on room temperature thermal coeffi- 
cients. Petit and Waddoups [11] developed a method of nonlinear analy- 
sis to failure. Symmetric laminates subject to biaxial loads were 
modeled. Loads were applied incrementally, and nonlinearities were 
introduced using stress-strain curves obtained from uniaxial tests. 

Good agreement was shown between theory and experiment for limited 
load conditions. Hashin, Bagchl and Rosen [12] and later Hashin, 

Rosen, and Pipes [13] used Ramberg-Osgood representations of unidir- 
ectional stress-strain data to model the nonlinear response of 
laminates. The deformation theory of plasticity was used, and inter- 
action of stresses for multiaxial stress states was accomplished via 
arguments based on isotropic J 2 theory [13]. 

Sandhu [14] used a method of fitting nonlinear lamina stress-strain 
data by cubic spline functions to predict nonlinear laminate response to 
failure. Loads were applied incrementally, and a flow rule of sorts was 
used to model stress interactions by introduction of an equivalent strain 
function. Compliances for successive steps were assumed to be functions 
of total equivalent strains. 

Hahn and Pagano [15] developed a method of including what might be 
referred to as second order thermal effects. They included the effects 


of stress- temperature coupling using an Instantaneous coupling coefficient. 
The consequences of omitting this effect were shown to be significant. 

In his modified lamination theory, Pagano [6,7] has extended 
Reissner's variational principle to laminated bodies. Using an assump- 
tion of stresses varying linearly with the through-thickness coordinate, 
he developed a set of differential equations along with the appropriate 
boundary conditions. His results compare well with the generalized 
plane strain finite element results of Wang and Crossman [16], Linear 
elastic material properties were assumed, and no thermal effects were 
included. 

More recently, Hashin, Rosen, and Pipes [13] developed a method of 
nonlinear laminate analysis including thermal effects. Simultaneous 
heating and loading were considered. As in the previous work of 
Hashin, et al_ [12], deformation theory was used. The results are 
startling in that they report experimental data that show that in the 
case of [0/90],. Graphite/Epoxy laminates the thermal expansion coef- 
ficients l ■* affected by time dependent (creep) characteristics and 
response. Results such as these appear to parallel works such as 
Mauri, et al_ [17], who report that moisture expansion coefficients are 
not only a function of temperature, but of the temperature at which 
the moisture was absorbed. 

The solutions referred to here as full field are so called because 
the problem is not limited by assumptions such as those for laminate 
analysis. Thus, the full set of stresses may be solved for, as well 
as the displacements for the complete mathematical model. These 
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solutions have been dominated by finite element work, but some early 
finite difference models are noteworthy. 

Pipes and Pagano [18] and Pipes [19] utilized a finite difference 
technique to solve the quasi-three dimensional elasticity equations 
for laminates. Interlaminar stresses were evaluated, and the capability 
for thermal effects was Included. Two mechanisms of interlaminar load 
transfer were discussed. More recently, Hsu and Herakovich [8] 
reported a perturbation solution to the quasi-three dimensional 
elasticity problem. The perturbation solution permits a better 
treatment near the "singularity" at the free edge than was previously 
available. The results of Hsu and Herakovich [8] are based on a 
stretching transformation of the boundary layer and are therefore 
dependent on an undetermined parameter of the transformation. They 
are thus qualitative in nature, but are valuable In that they place 
oounds on the region of applicability of the standard finite 
difference solution, as well as showing relative shapes of stress 
profiles in the boundary layer. Thermal effects were not included. 

Isakson and Levy [20] used the finite element method (FEM) to 
model a laminate in a state of plane stress. The plies were modeled 
as orthotropic layers separated by isotropic shear layers. Inter- 
laminar normal stress was omitted. The technique of using isotropic 
shear layers had been previously used by Puppo and Evensen [21], but 
not in conjunction with the FEM. Subsequently, Levy, Armen, and 
Whiteside [22] extended the work of Isakson and Levy for plastic 
deformation of the shear layer. A Ramberg-Osgood representation for 
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boron/epoxy was used as a model for the nonlinearity of the shear 
layer. Results were presented for Interlaminar stresses around holes, 
and the effects of stacking sequence were discussed. No thermal effects 
were considered. 

Herakovich and Brooks [23] used the FEM to examine stress distribu- 
tions In uniaxial and cross-piled laminates subjected to thermal and 
uniaxial strain loading. Material models were linear elastic. Inter- 
laminar stresses were shown to be significant near free edges. 

Herakovich [24] later showed the effects of stacking sequence on inter- 
laminar thermal stress. His analysis Indicated an effect on the overall 
strength of the laminate. The analysis was again linear elastic. Lin 
[25], and later Dana [26], and Dana and Barker [27] used a 72 degree 
of freedom, three dimensional orthotropic isoparametric finite element 
element to model laminates. All results were linear elastic, and 
thermal effects were not included. Stresses were computed near holes 
and free edges. Comparisons of hole shape and size were made with 
regard to Interlaminar stresses. 

Rybicki has used several variants of the FEM to examine the 
behavior of composites. Complementary energy methods [28], stress 
function methods [29], and standard displacement techniques [30,31] 
have been used. More recently, Rybicki and Schmuesser [32] used the 
three dimensional finite element program SAP IV to predict stress 
distributions near holes in linear elastic laminates. Stanton has used 
a parametric cubic modeling system [33,34] to analyze nonlinear material 
behavior in laminates. His analysis is three dimensional, its basic 
drawback being that the complexity of the element and resulting large 
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bandwidth dictate the use of a minimization technique for solution. 

Foye [35] has used a micromechanics approach based on the FEM to 
model nonlinear behavior of composites. Including both nonlinear time 
Independent and time dependent effects. By using a micromechanics 
approach and thus Introducing fiber and matrix properties Independently, 
the difficulty of orthotropic plasticity as noted by Hill [36] Is 
avoided. Effects of differences In cure cycle are examined. The cure 
cycle for these analyses was approximated by a series of finite steps 
of temperature followed by a constant temperature until the next step. 

Plfko, Levine and Armen [37] have developed a method of three 
dimensional stress analysis of composites. Their work is based on the 
FEM displacement formulation. Plasticity Is Included via the normal 
flow rule (Drucker's postulate) and an incremental analysis. Hardening 
is either nonlinear or linear, with nonlinear hardening introduced via 
Ramberg-Osgood approximations. Kinematic hardening [38] is used. 
Unfortunately, the kinematic hardening coefficients are based on 
heuristic arguments, and are not fully checked out [39]. 

Renieri and Herakovich [40] used a quasi -three dimensional finite 
element analysis to model the response of laminates to thermal and 
mechanical laodlng. Nonlinear material properties were introduced 
via one-dimensional Ramberg-Osgood representations. Interaction of 
stresses was not considered, and loading and unloading were along the 
same path. Hence the analysis may be called nonlinear elastic. 

Differing properties In tension and compression were included. 

Various non-displacement finite element formulations have also 
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been used for stress analysis of composites. Lee [41] and Wang, et al_ 
[42,43] have used the hybrid stress formulation of Flan [44] to 
perform finite element analyses of composites. Several problems exist 
with this approach. In the boundary layer, steep gradients of stress 
are known to exist. Thus, In the hybrid stress approach one must 
Include significant numbers of higher order terms In the assumed stress 
polynomials. The analysis Is thereby complicated by a large number of 
unknowns per element. Kathlresan [45] and Atlurl and Kathlresan [46] 
have used the hybrid displacement formulation for modeling of composites. 
Such an approach appears to have favorable qualities for Investigation 
of fracture mechanics and laminates containing Initial flaws, although 
treatment of nonl Inearl ties may be difficult. 

In summary, considerable work has been done In lamination theory 
studies and finite element studies. In lamination theory, the work of 
Pagano [6,7] stands out as the only one capable of Inclusion of Inter- 
laminar effects. Problems will undoubtedly arise, however, If 
extension to nonlinear material behavior and/or thermal effects is 
attempted. The finite element method has been used in all forms to 
solve laminate problems. The current state-of-the-art appears to be 
two dimensional (quasi-three dimensional) analysis with nonlinear 
elastic material properties (different In tension and compression) 
with first order thermal effects or three dimensional analysis with 
nonlinear material behavior. The nonlinear three-dimensional analysis 
of Stanton [33,34] apparently has not been widely accepted by the 
composites community. The areas which seem to be lacking are two 
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dimensional (quasi-three dimensional) analysis with time dependent 
material properties and three dimensional analysis with a full 
plasticity treatment and moisture and creep effects Included. No 
finite element analysis to date appears to have Included second order 
thermal effects as described by Hahn and Pagano [15], even though the 
formulations have been made by Ueda and Yamakawa [47]. 

When considering the modeling of a material system, one must 
always survey the availability of material property data. For 
composites, the data package Is far from complete. Quasi-static stress- 
strain data at room temperature Is relatively abundant. The most 
effective means of Implementing temperature dependence seems to be the 
percent retention at temperature curves as used by Renlerl and 
Herakovlch [40], even though their validity is only for certain 
strain ranges. Time dependent properties such as viscoelastic moduli 
and creep compliances are only now being generated [48] for even the 
most common composites. Temperature dependence of these properties Is 
only known as far as tests at a very limited number of temperatures. 
Effects such as the stress dependence of thermal expansion coefficient 
noted by Hashln, et al^ [13] and the absorption temperature dependence 
of moisture coefficients noted by Mauri, et al [17] are only pointed 
out as anomalies. These effects are likely the result of the complex 
polymerization characteristic of resin matrix composites. Understanding, 
and ultimately Intelligent and accurate modeling, must await more 
knowledge of the physical processes of the materials Involved. 


Chapter 3 

NONLINEAR FINITE ELEMENT FORMULATION 


3.1 Introduction 

The finite element used for the nonlinear studies Is the 24 node 
Isoparametric element described In Appendix A. The basis for an 
Incremental solution of the nonlinear load-deflection equations Is 
also presented In Appendix A. This chapter details the pseudoloads 
used for Implementation of plasticity using the method of Initial 
strains, as well as the nonlinear thermal and plasticity modifications 
to the elastic constitutive relationships. It should be emphasized 
that the thermal and moisture effects are uncoupled . That Is, the 
temperature and moisture distributions are known a priori by means of 
some also uncoupled thermal or moisture diffusion analysis. Further- 
more, the plasticity developed here Is time Independent. No visco- 
elastic or creep effects are considered. However, using the Initial 
strain method for time dependent effects proposed by Zlenklewlcz and 
Cormeau [49], their Inclusion Is straightforward. Several additional 
assumptions are necessary In the developments to follow. These will 
be stated as they are Introduced. 

Several types of nonl Inearl ties exist for laminated composites. 
The constitutive parameters are temperature and moisture depen- 
dent, history dependent (plasticity) and time dependent (creep). 

The thermal expansion coefficients are temperature and moisture 
dependent. Moisture expansion coefficients are temperature dependent. 
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Various "coupling" of these effects may be shown to exist. Hahn and 
Pagano [15] have used a lamination theory model to show the significance 
of the so-called "coupling" terms which exist between temperature ef- 
I fects and stress. In the present work, these have been termed "second 

order" thermal effects and are Included In the constitutive relation- 
ships. 

In the developments of this chapter, contracted tensor notation for 

i | 

, strain (e^) will Imply engineering shear strains, while full strain 

i tensor notation < e 1j' Implies tensor shear strain. Conversion from 

tensor to engineering shear strain Is assumed to have been performed 
i as required. 

I 

> 

3.2 Incremental Solution of Nonlinear Equations 

The method chosen for the Incremental solution of the nonlinear 
load-deflection equations is referred to by Desal and Abel [50] as the 
direct Incremental Initial strain method. All nonlinearities are 
Introduced as pseudoloads, with corresponding Initial strain modifica- 
tions to the constitutive relationships. One of the benefits of the 
Initial strain method is that the global stiffness matrix need only 
be assembled once. Furthermore, unless a change of specified 
displacement boundary conditions occurs, the solutloi. requires only 
one triangularlzatlon. The utility and computational efficiency of 
the initial strain method Is widely recognized. Pifko, et a| [37] 
used an Initial strain method in the development of the PLANS system. 

The following paragraphs contain a description of the initial strain 
method as Implemented for the present study. 

I 
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3.2.1 Direct Incremental Initial Strain Method 
The starting point for the method Is the Incremental load-deflection 
solution of Appendix A, 


K ij (Aq j )k+1 " < A( M k+1 * < A0 ? J ) k+1 ♦ UQ? ,P ) k ♦ (3-D 


where 

k 

k+1 

Uq /* 1 

(AQ 1 ) k+1 

UQ ° J ) k+1 

UQ?* p ) k 

(R 1 ) k 

N 


1 • 1,2 N 


* previous load Increment 

■ curr*nt load Increment 

z current Incremental nodal displacement vector 

* mechanical load Increment 

■ Incremental pseudoload due to temperature effects 

■ Incremental pseudoload due to plastic flow during 
previous load step 

* residual from total equilibrium at end of previous 
load seep 

* global elastic stiffness matrix. 

* total number of degrees of freedom. 


Note that the method Is directly applicable to linear elastic problems, 
where only one load str:* to full load results In a purely linear 
analysis, using the linear elastic constitutive relationships. 

The basic procedure Is a direct stepping technique, where, given 
the state variables and material properties at the beginning of a step, 
the Increments are found by solution of Equ's. (3.1) and the appro- 
priate constitutive relationships. Total quantities are computed by 
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addition of Increments to previous totals. This procedure Is repeated 
until the full load Is applied. 

3.2.2 Computation of Incremental Initial Strain Pseudoloads 

For the Initial strain effects, the Incremental pseudoload Is 

given by 

aqJ • ffl B^D^AeJdn (3.2) 

where Bjj, D^j, and n are the same as described In Appendix A, and Ac® 

Is the Increment of initial strain. For plasticity, Ac® Is the 
plastic strain Increment for the previous load step. For temperature 
effects, Ac® Is the Increment of thermal expansion for the current load 
step. The computation of these terms will be discussed later. The 
volume Integral of (3.2) Is computed by Gaussian quadrature, as 
described In Appendix A. The elasticity matrix Is evaluated at 
the same temperature as the qlobal stiffness matrix. 

3.2.3 Additional Assumptions 

Several additional assumptions are made for implementation of the 
direct Incremental Initial strain method. All strains and displacements 
are assumed to be small. Furthermore, the anisotropy (transverse 
Isotropy for the present study) ;s assumed constant for the entire 
analysis. That Is, angular orientations do not change during the 
analysis, nor Is the anisotropy altered by plastic flow. As previously 
stated, all thermo-elastic, thermo-plastic, hygrothermal , etc., effects 
are actually uncoupled, the present analysis computing the elasto- 
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plastic effects of thermal and r.ioisture fields known a priori from 
some other analysis, also uncoupled. It is also assumed that the 
total strain can be written as the sum 

e i * e i + e 1 + £ i + £ i * * “ (3.3) 

where contracted tensor notation is used and e^, i = 4,5,6, are 
engineering shear strains, and 
eT = total strain 

A 

e.. = elastic (or mechanical) strain 

a 

e.j = thermal expansion strain 
M 

e" = moisture expansion strain 

P 

e.j = plastic strain. 

This assumption is standard in continuum analysis, and is necessary to 
facilitate the computation of the mechanical strain. 

3.3 Incremental Nonlinear Constitutive Relationships 

The derivation of nonlinear constitutive relationships which 
accurately represent the physical behavior of laminated composites and 
which are consistent with pseudoloads to be applied to a finite element 
model is the crux of the present study. Efforts have been made to 
utilize parameters which are common within composites terminology, such 
as the Ramberg-Osgood technique for representation of nonlinear stress- 
strain relationships. Certain of the assumptions made in the present 
study are commonly used. The method of including second order thermal 
effects in the finite element context is due to Ueda and Yamakawa [47], 
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with orthotropy modifications as suggested by Yamada [51]. Second 
order thermal effects, although rarely considered, were shown to be 
significant by Hahn * . I Pagano [15]. 

3.3.1 Nonlinear Thermal Effects 

As previously stated, thermal nonlinearities are cf two types. 

The instantaneous thermal expansion coefficient, o^, defined by 

3e 1 

“1 = “3T (3>4) 

Is dependent on temperature. Also, the "elastic" coefficients of the 
material are temperature dependent. In the present study, the first 
effect has been termed "first order", while the effects of the latter 
on thermal stress have been termed "second order". The following 
paragraphs describe the first and second order thermal effects and 
the corresponding elastic constitutive relationships. 

3. 3. 1.1 First Order Thermal Effects 
The Incremental form of Equ's. (3.3) is 

ArJ = At® + At® + Ac^ + Ac^. (3.5) 

In the absence of moisture and plastic flow, Equ's. (3.5) (after 
rearranging) become 

At* » Atj - Ar®. (3.6) 

For first order nonlinear thermal effects, the incremental thermal 
strain is approximated by 
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Ac® - c^aT (3.7) 

where AT is the temperature change over the increment and is an 
average coefficient of thermal expansion over aT. Note that if the 
a,j are not functions of temperature, Equ's. (3.7) reduce to the familiar 
constant property relationship of Appendix A. The Incremental constitu- 
tive relationships are 

4 °i = (3 - 8 > 

where Ao^ is the stress increment, and is the elasticity matrix 
(rotated to account for orthotrooy if required) for the material, and 
is evaluated at the mean temperature over aT. Combining (3.6), (3.7), 
and (3.8), we obtain the constitutive relationship for first order 
thermal effects: 

Ao i = D ij^ Ac j ' a j AT ^ (3.9) 

The first order incremental thermal pseudoload, aq9, is, using Equ's. 
(3.2) and (3.7), 

AQ® = fff B j 1 D j k a k AT dn. (3.10) 

Equ's. (3.9) and (3.10), along with the strain-displacement relationships 

At i = (3.11) 

form the basis for the first order nonlinear thermal analysis of the 
present study. 
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3. 3. 1.2 Second Order Thermal Effects 

An effect shown by Hahn and Pagano [15] to be significant for 
composites, where material properties may change rapidly with tempera- 
ture, and by Ueda and Yamakawa [47] for welding problems where large 
temperature gradients exist, Is the "coupling" between total stress 
and thermal stress. Note that the analysis Is not truly coupled thermo- 
structural, but that mutual influence of effects Is being considered. 
Suppose that the mechanical strain, e®. Is a function of temperature 
and total stress, or 


1 1 = 


(3.12) 


From the constitutive relationships. 


e _ c 

■I - s ij”o 


(3.13) 


where S^ is the elastic compliance matrix, or 


S ij - ( D ij> 


-1 


(3.14) 


Now, realizing that the are functions of temperature, the differen- 
tial elastic strain becomes, from Equ's. (3.13) 


df i = s ij d °j * TT* °0 dr> 


(3.15) 


or. In Incremental form 


' f 1 = 


aSi . 

S 1j A "j + Tt °j AT - 


(3.16) 
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The se:ond term in (3.16) is what Hahn and Pagano [15] refer to as the 

"coupling" term. Note that S,,. is evaluated at the mean temperature 
3S.. J 

over AT, as is . The incremental constitutive relationship becomes 


ta i * V E j - D ij 



(3.17) 


Note that in Equ's. (3.17) the incremental stress Ao^ is a function of 
total stress a k « For accuracy, a k should be evaluated at the midpoint 
of the load (temperature) step, but it is only known at the beginning 
of the step. The incremental solution thus is implicit, and an 
iteration or some additional approximation is required. In addition, 
the initial strain necessary for computation of the incremental thermal 
pseudoload is 



and the incremental pseudoload becomes 

j ATdrt. (3.19) 

Equ's. (3.17) and (3.19) are the constitutive relationship and thermal 
pseudoload, respectively, for the second order nonlinear thermal 
analysis. 

3. 3. 1.3 Comments on Thermal Nonlinear Analysis 
Several assumptions have been made in the thermal nonlinear 
formulations. They include the assumption of constant instantaneous 
thermal expansion coefficients and changes of with temperature over 


«?■ 


Iff 


B . .D.i 
Ji Jk 




3S 


k£ 

3T 
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an entire temperature step. The accuracy of these assumptions in 
practice Is dependent on the rates of change of the coefficients of 
concern and the size of the temperature step. Where steep gradients 
exist, small temperature steps may be required, as opposed to large 
steps being acceptable where property behavior is less temperature 
dependent. When including second order nonlinear thermal effects, the 
magnitude of total stress is also of concern. As with all nonlinear 
solutions, considerable care must be taken In choosing load steps, or 
the software used must do considerable checking to assure that the 
assumptions are not violated. Such checking Is computationally 
expensive and seems to be done rarely If at all. 

3.3.2 Plasticity 

The plasticity used In the present study is based on the yield 
function, equivalent stress, and equivalent plastic strain increment 
derived for anisotropic materials by Hill [36]. It should be noted 
that for isotropic materials Hill's relationships reduce to the familiar 
von Mises (distortion energy) theory. Inherent in this formulation 
is the incompressibility of plastic strains and the assumption that 
hydrostatic (spherical) stress states result in no yielding or plastic 
deformation. It has been noted [52] that this is not the case for all 
composite materials. Hill's yield function, or some variant, is 
nevertheless widely used. Once yielding is detected, a mathematical 
relationship for progression of plastic flow is required. The most 
commonly used relationship is Drucker's postulate [53], which states: 
"Suppose an external agency to produce a non-zero displacement by 
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adding a set of loads to a time- Independent system In static equilibrium. 
The system Is stable If the work done by the external agency Is 
positive for all permissible added loads." Mathematically, Drucker's 
postulate Is common ,y combined with the so-called "normality" princi- 
ple and written as 


dc 


P 

ij 


3 °ij 


(3.20) 


where 

p 

dc^j = plastic strain Increment 

f( 0lj ) ■ yield function (associated plasticity) 

"u ■ stress 

x = nonnegative proportionality constant. 


When Drucker's postulate Is used In a numerical scheme, such as an 
Incremental solution, several assumptions are made. The numerical 
value of x Is assumed to be constant over the entire step. Equations 
(3.20) are often referred to as the "normal flow rule". If we consider 
f ( o . j ) to be a yield hypersurface In nine-dimensional stress space, then 

D 

Equ's. (3.20) require that the vector de.,< be "normal" to the hyper- 

1J 3f(o..) 

surface. The assumption Is thus also made that — — , or the dir- 

P ’"u 

ection of the de^ vector, is also constant over the increment. As 
with all incremental schemes, plasticity for the case of nonpropor- 
tional loading is at best convergent in the limiting case of 
infinitesimally small increments. 

It should be emphasized that Drucker's postulate is only a 
postulate, not derivable from first principles. It is a tool which 
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facilitates computation of plastic strains based on stress state. The 
proportionality constant, a, Is usually derived from plastic work 
considerations [36,38]. Other Investigators [54,55,56] have chosen 
approaches In which Drucker's postulate Is not used, often with good 
results. Drucker's postulate is used In the present study, with the 
full realization that it is purely an assumption, and that other 
means could be utilized. 

3.3.2. 1 Incremental Plastic Pseudoload 

The incremental pseudoload for use In the (k + l) th step is com- 

P th 

puted using the plastic strain increment, Ae^ derived from the k 
step. Using the standard pseudoload form of Equ's. (3.2) in the nota- 
tion of Equ's. (3.1 ), 

(AQ?’ P ) k+1 = /// B j1 D jl (Ac P ) k d». (3.21) 

The volume integral is omputed, as before, using Gaussian quadrature. 
Thus, the incremental pseudoload is computed using the prior step 
plastic strain increment, which, as shall be shown, is derived using 
some values from the (k - l) th increment. The scheme is seen to be 
increasingly dependent on the choice of load increment, the smaller 
the load increment the better until round-off and discretization error 
become of the same order as the error due to linearization. 

3. 3. 2. 2 Nonlinear Thermoplastic Constitutive Relationships 

Once the total, plastic, and thermal strain increments are computed 
for a given load step, c constitutive relationship is required to compute 
the stress increment. The following developments follow those in the 
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text by Hill [36], and are restated here for completeness. It should 
be noted that all stresses, strains, yield conditions, etc., must be 
rotated to the principal material axes. The notation adopted Is s^ 
and e^ for stress and strain, respectively, In material coordinates. 

The 1-dlrectlon will be assumed to be parallel to the fibers, the 
2-dlrectlon to be transverse, In the plane of the thin laminae, and 
the 3-dlrectlon normal to the plane of the lamina, as shown In 
Figure B.3. 

The yield function of Hill [36] for Initial yield Is 

f (s 1 ) * F(s 2 - s 3 ) 2 + G(s 3 - Sj) 2 + H(Sj - s 2 ) 2 + 2 Ls 23 

+ 2Ms^ 3 + 2Ns 2 2 * 1 , (3.22) 

where 


7* 

1 

?' 

i 

X 2 

(3.23a) 

7 + 

i 

z 2 ' 

1 

V 

(3.23b) 

7* 

1 

V 2 ' 

i 

7 

(3.23c) 

2L = 

1 

7 


(3.23d) 

2M = 

i 

7 


(3.23e) 

2N = 

i 

7 


(3.23f ) 
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and 

X, Y, Z, R, S and T are yield strengths In the principal material 
coordinate system directions 1, 2, 3, 23, 13, and 12, respectively. 

Equ. (3.22) Is valid, for a hardening material, only for Initial yield. 

As yielding progresses, so does the numerical value of f(s^). 

The model used In the present study Is known as Isotropic hardening. 
The basic premise of Isotropic hardening is that as the yield hyper- 
surface grows In stress space. It grows uniformly in all directions. 

As an example, consider the von Mlses ellipse In two dimensional stress 
space for an Initially Isotropic material (Figure 3.1). The solid 
ellipse represents f(o^) for Initial yield, and the dashed ellipse 
represents f(o^) after some plastic flow has occurred. In either case, 
a stress state In or on the ellipse is elastic, while stress states 
outside the ellipse are governed by the normal flow rule. Cyclic 
loading might result in loading into the plastic range, followed 
by elastic unloading and reloading to subsequent yield followed by 
additional plastic deformation. 

Isotropic hardening has two serious shortcomings. The plastic flow 
governed by the normal flow rule expands uniformly (or isotropically) in 
stress space. Thus, an initially isotropic material remains isotropic, 
and an initially anisotropic material undergoes no change in its 
anisotropy. Experimentally, this has been shown to be inaccurate. An 
example is cold rolled sheet which is initially isotropic, but exhibits 
anisotropy following the plastic flow of the rolling process. The other 
deficiency of isotropic hardening is that it allows for no Bauschinger 
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FIGURE 3.1 VON MISES ELLIPSE FOR INITIAL AND SUBSEQUENT YIELD 


27 


effect, a commonly observed phenomenon. When using the flow rule, the 
alternatives to isotropic hardening are kinematic hardening [37], In 
which the yield surface translates, but does not grow. In stress space, 
or some mixed model [57] which Is some mixture of Isotropic and kinematic 
hardening. Both models exhibit computational difficulties, especially 
for anisotropic materials. 

Also used In the elastoplastlc analysis Is the equivalent stress, 
as given by Hill [36], 

* * 2(F+6+H) " s 3^ + G( s 3 - *,) + H(s^ - Sg) 

+2Ls| 3 + 2Ms* 3 + 2Ns^]. (3.24) 

Note the difference In the definition of s in Equ. (3.24) and f(s^) in 
Equ. (3.23). Realizing the basic similarity In the definitions, it is 
assumed that yielding begins when 

5( Sl ) * s(e P ,T) (3.25) 

where the functional notation indicates that s Is a function of the 
current stress state, and s is a function of the equivalent plastic 
strain (history) and temperature. Taking the differentials of both 
sides, and realizing that 5 and s have the same functional form, 

^ dSlj - -jfr di” ♦ ff dT (3.26) 

where the equivalent plastic strain Is related to the equivalent plastic 
strain Increment as defined by Hill [36], 


J 
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de P • */| (F+G+H) [ 

4 H(Fde^-Gde22) 
(FG+GH+HF) 2 


FfGde^-Hde^) 2 ♦ GCNde^-Fde^ ) 2 
(fg+ghVhf) 2 

(de 2 P 3 ) 2 (dt P 3 ) 2 (de P 2 ) 2 I 

+ "TC — + “TfH — “”7R — J 


(3.27) 


where e 


P 

1J 


. 1 M. are tensor shear strains, or 


e 


1j 


1 

7 


3Uj 3U. 

( — - ♦ — 1 “•) 
' 3Xj 3X^ ' 


1 ,j»l ,2, 3, 4, 5, 6. 


(3.28) 


where and are the standard displacement and position vectors 
of continuum mechanics, respectively. 

It should be noted that Hill's equivalent plastic strain Increment 
(Equ. 3.27) Is, In general, only defined as an Increment. It cannot be 
directly Integrated unless some assumption Is made of the loading 
history. The often used deformation theory of plasticity utilizes the 
assumption of proportional loading, in which loads and stresses Increase 
In proportion, to facilitate Integration of the equivalent plastic 
strain Increment and the flow rule over the entire load history. The 
assumption of proportional loading Is unacceptable for cyclic loading. 
Indeed, for a composite material cured at an elevated temperature, 
cooled to room temperature (resulting In some compressive residual 
stress), and then loaded In tension, any load history Is cyclic to 
some extent. It therefore appears that deformation theory Is unaccep- 
table, unless applied In some quasl-lncremental fashion, to problems 
Involving curing and subsequent loading of composites . The incremental 
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or flow theory of plasticity will be used In the present study. 

Several comments are necessary concerning Equ. (3.26). The func- 
tional form of the effective stress, 5, Is given by Equ. (3.24). The 

terms -r—- are computed using the closed-form partial derivative. The 
3S 1 

anisotropy coeffclents F, G, H, l, M, and N are functions of temperature 

which may be determined experimentally, although not without difficulty. 

Thus |y reduces to terms of the form |y, etc. The terms remaining are 

and de P . The term de P Is Hill's equivalent plastic strain increment, 
3 * a? 

Equ. (3.27). The term — will be referred to as the hardening coeffl- 

cient. Hill [36] and later Plfko, et a1_ [37] have conmented on the 

difficulty of Its determination. Mathematically, -t- may be thought of 

-P 3 * 

as the slope of a plot of s vs e . Unfortunately, for anisotropic 
materials it Is not uniquely determined for a given stress-strain state, 
but Is a function of the complete history of the material. 

The first term of Equ. (3.26) determined Is the equivalent plastic 

-P P 

strain Increment, de . The increment of plastic work, dW , done by a 

p 

plastic strain increment de^ at stress s^j is 

dW P * S^deJj. (3.29) 

Similarly, as shown by Hill [36], 

dW P * sde P . (3.30) 


Multiplying both sides of the normal flow rule (Equ's. 3.20) by s^ 

o 3f(s<J 


s 1j de ij ’ s ij X “as 


ij 


(3.31) 
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Using (3.29) and (3.30) in (3.31), 


3f (S 1 . ) 

5 ^ = V “^ij" 




or 


6 f,,lAl!pA 

5 9s ij 


Using Equ. (3.3; in Tqu. (3.26), 


(3.32) 


(3.33) 


jJL ds = _ii S J1 + a£ d T 

3s ij d# « 3i P 5 3s ij 3T 


(3.34) 


Writing the separation of total strain (in material coordinates) into 
components of Equ's. (3.3) in incremental form results in 

Ael = Ae® + Ae? + Ae P , (3.35) 

where all strains are engineering strains. 

Realizing that 

As. = Ae®, (3.36) 

where D?. is the elasticity matrix, Equ's. (3.35) result in the incre- 

' J 

mental constitutive relationship 

< 3 - 37 > 

where Ae® is given by Equ's. (3.18) as 

Ae i = (a i + 


(3.38) 
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where S^j is the compliance matrix, is the current stress, T the 
temperature, and the instantaneous thermal expansion coefficient. 

p 

The plastic strain increment, Ae;., is given by the flow rule, Equ's. 
(3.20) to be, in incremental form, 

3f(s.) 


Ae i = x as. 


(3.39) 


Writing Equ's. (3.34) in incremental form and noting conversion to 
engineering shear strains, 


as _ , al s i as T 

“T 1 ' w 


(3.40) 


Substituting (3.37), (3.38), and (3.39) into (3.40) 


aS 

as n e _ as n e / x jk . as 


da n c da n c / , 

5sT - sij °ij ( “j * 


3T s k' 4T ' ft 4T 


ji n e . . . _js !i 

35, 1j 3 Sj si P g 35, 


(3.41) 


where Ae.. has been used to represent the total strain increment. 
Solving (3.41) for x. 


-iS-D.-Ae. - -5 s 
3s, "ij“j as, 




as at 

5T 4T 


35 S i 3f < s i> 
si* 5 s 


as j 


+ -ilD*. 

as i ij 


af(s.) 

~*T 


(3.42) 


where and S^ are evaluated at the current temperature. Note that 
aside from being the proportionality constant in the flow rule, X serves 
as a loading indicator. The three conditions are 
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X < 0 

(3.42a) 

X ■ 0 

(3.42b) 

X > 0, 

(3.42c) 


which are denoted a; unloading, neutral loading, and loading, respectively. 
Unloading (x<0) Implies that no plastic flow is occurring, the stress 
state being inside the yield hypersurface, and the elastic constitutive 
relations are used. Neutral loading (x*0) implies that no additional 
plastic flow is occurring, but the stress state is moving along the 
yield hypersurface, and the elastic constitutive relations are used. 

Loading (x>0) implies that additional plastic flow is occurring, and 
the yield hypersurface Is growing in stress space. The plastic 
constitutive relations are used. 

In the present study, the numerical value of Hill's yield condition 
(Equ. 3.22) is computed at the end of each load step for each Gauss 
point. If unloading Is detected, the elastic constitutive relations 
are used until the maximum value of Hill's yield condition for that 
Gauss point is reattained. Thus, cyclic loads may be treated. 

The Incremental plastic constitutive relationship becomes 

e e , aS ik v e 3f(s i J 

«i ' D V e j - * "# s k> 4T - XD ?j -sr- < 3 -«> 

J 

where x is given by Equ. (3.42). In the case where a transition from 
elastic to plastic is made during a load step, some adjustment to 
Equ's. (3.43) is required. This region has been termed the transition 
region. Marcal [58] gives a computationally efficient procedure for 
its treatment. Suppose that a fraction, denoted by m, of the load 
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step was elastic. Then the substitution of 

X = (1 - m)x (3.44) 

Into Equ's. (3.43) reduces the plastic strain, and thus corrects the 

transition region stress increment. 

3. 3. 2. 3 Nonlinear Orthotropic Hardening Coefficient 

The remaining term occurring in the plasticity developments is the 

hardening coefficient, -^t-. As previously stated, this quantity has 

ae 

been pointed out as being difficult to determine for anisotropic 
materials. Biffle [59] used an assumption of equal plastic work in 
all directions. Pifko, et al [37] used heuristic arguments based on 
the method of weighted averages. Such approximations are necessary, 
but are obviously postulates. The technique developed for use in the 
present study utilizes assumptions which are also heuristic but will be 
shown to yield reasonably accurate results for both isotropic and 
orthotropic materials. 

Suppose, in an incremental analysis, that the state variables are 
known at the end of an increment designated as k. For initial yielding, 
the step k will have been elastic. For continuing plasticity, the state 
will have been computed using the plastic constitutive relationships. 
Also suppose that the uniaxial stress-strain relations, in principal 
material directions for orthotropic media, are given in Ramberg-Osgood 
[60] form, 

e i = fT + 6 i (s i) ni i=1 ’ 6 (3.45) 
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where the summation convention has been suspended for the entireity of 
this development, is Young's modulus for the i direction, and 8^ and 
n. are known as the Ramberg-Osgood coefficients for the i direction. Let 
the stress increment for the k step be represented by (as^) . We 
define the max norm [61] of the stress increment as 


I | ( as . ) k | I = max | (as . ) k | 
1 1 < i <6 1 


(3.46) 


The normalized stress increment, Us^) k , is defined by 




(as/ 

l|(as) k || ' 


(3.47) 


The model developed in this study for computation of -^4 utilizes 
the normalized stress increment for the k cn step and the total stresses 
at the end of the k^ step. For the purposes of computing the hardening 
coefficient only, proportional loading is assumed during the (k+l)** 1 
step. Hashin, et a[ [12] have used proportional loading for the entire 
loading history, and have discussed its accuracy for loading paths 
which are "neighbors" of proportional loading. Hence, the assumption 
of proportionality from step to step is reasonable unless steps are 
large or contain load reversals. Note that step-to-step proportionality 
is used only for computation of the hardening coefficient. Once it is 
computed, the previously developed plastic constitutive relations are 
used to compute the stress state at the end of the (k+l)^* 1 step. 


As done by previous investigators [37,40], the nonlinear, or plastic. 
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uniaxial strain Is assumed to be the power law portion of the Ramberg- 
Osgood approximation, or 

Oj 


Differentiating, 


P 1 

e,j - B i (s i ) . 


P 

de. = B i n i (s 1 ) ds.., 


(3.48) 


(3.49) 


where the assumption Is made that plastic strains are related to 
changes in stress. Inherent in the use of a Hill-type criterion is 
the incompressibility of plastic strains, or 


i=l 


e 1 = 


(3.50) 


This condition must be strictly enforced for consistency. 

Consider the case of a prismatical bar of an isotropic elastic- 
plastic material subjected to an axial load. When loaded in tension 
into the plastic region, the axial plastic strain is related to the 
stress through some relation such as Equ. 3.49, where e and n are the 
appropriate Ramberg-Osgood parameters. Transverse plastic strains 
occur such that the incompressibility condition (3.50) is satisfied, 
but stress changes are small, such that (3.49) are not satisfied in 
the transverse directions. We therefore speak of the axial direction 
as being the "dominant" direction of plastic flow, and use the axial 

p 

Ramberg-Osgood parameters. In this manner, we obtain d e ax ^ a i* From 
the flow rule. 
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de 


1j 


; af(s 1J } 


as 


1J 


(3.51) 


p p p p 

Indicating the existence of the ratios de£:dej and de 3 :de^. Therefore, 

p p 3f(s^i) p 

de 2 and de 3 may be computed given 3 — ■ and de^ Note that no such 

conditions exist on plastic shearing strain. It Is therefore assumed 

that the shear stress-shear strain relations follow the Ramberg- 

Osgood curves determined for pure shear. 

Now consider the generalization of the previous arguments to an 

orthotropic material. The "dominant" direction Is determined as the 

3f(s^) 

direction corresponding to the maximum absolute value of 


DS 


ILL 


ij 


Ramberg-Osgood coefficients govern the flow In the dominant direction, 

3S 


and transverse flow is determined by the relative values in 


as 


and 


1j 


the incompressibility condition. Plastic shearing strain follows the 
Ramberg-Osgood approximations determined by pure shear tests. 

For purposes of illustration, consider the case where the 
1-directlon has been determined to be dominant. Then 


and 


P n l 

de 1 * B^Sj ds 1 


P n l 

de 2 = c 2 e 1 n 1 s 1 ds 3 


(3.52) 


(3.53a) 


n r l 


de 3 • CjB,n, Sl ds, 
3f(s,,) 

where c 2 and c 3 are determined from — — . 

1 *1 


(3.53b) 


Furthermore, assume that 


the stress increment, stress, and plastic strain increment are known 
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at the end of the step, and that the loading Is proportional from 
the k th step to the k+1^ step. The projected stress Increment for the 
k+l th step Is thus 

(As) k+1 « * k (As) k (3.54) 

where $ Is the proportionality constant for the step k to k+1. Sine#* 
the effective stress Is a function of the Individual stresses, the 
differential change in the effective stress is 


d5 ■ aif ds l + jq ds 2 + ds 3 * 9 I 7 ds 4 + ds 5 * af| d V 


(3.55) 


Since the stresses and Increments for the k th step are known, using 
the proportionality assumption, the increment of the effective stress 
for the k+l th step may be estimated as 


(as) 


k+1 - A k x Ac . as . . 3i Ac . as . a 5 Ac , as Ac ,k 

" ♦ { a?7 AS 1 + aij As 2 + ai^ as 3 + al^ as 4 + a?; As 5 + al^ As 6 } 


(3.56) 


where Incremental changes have been substituted for infinitesimals. 

The next operation is to project the effective plastic strain incre- 
ment of Equ. (3.27) from known quantities for the k^ step to the k+1*^ 1 
step. Using (3.52), (3.53), the differential shear stress-shear strain 
relations similar to (3.52), and the proportionality assumption, the 
effective plastic strain increment becomes 


P k+1 k /p P] - ! k ? F(GCp-HC^) 

(de P ) k+1 - (F+G+H) {[B-n-s, 1 (As,) k ] Z [ *— *- 

J 111 1 (FG+GH+HF) 
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G(HC 3 -F 2 +H(F-GC 2 ) 2 n n 4‘\ A . %k,2 (3 57 \ 

♦ 2 5 -*—] + gr [e 4 V4 (a V ] 1 } 

(FG+GH+HF) 2 BL 

+ SR te 5 n 5 s 5 5 + 5R [e 6 n 6 s 6 6 

where the difference In tensor and engineering shear strains has been 
accounted for. Using the normalized stress Increment of Equ. 3.47 and 
realizing that the hardening coefficient (-Jp) may be approximated 
by (-^|) k+1 . we find that the proportionality constant and the max 


of the prior step stress Increment cancel out, leaving 


- k+1 ,r k+1 

(-4) » (-rp) 

Ae 


where 


■ £ <“/ * 4 ( “2 )k * 4 ( “3 ,k 


* £ (“/ ♦ 4 l “s’ * 4 


, n,-l . , 7 F(GC 2 -HC 3 ) 

<« p ) k+1 4^ ( t e i n i s i <“# 3 


G(HC,-F 2 +H(F-GC 3 ) Z ! V 1 /.: \k-i 2 

♦ 5 X 2 ” 3 + 8U [ W4 AS 4 3 

( FG+GH+HF r 


+ gR [6 5 n 5 S 5 5 (AS 5 )k]2 * C W6 {AS 6 )k]2}? ’ 
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Note that In Equ's. (3.58) and (3.59) the stresses, s^ and effective 

stress gradients, r— are computed at the end of the k th load step. 

5 i 

Several observations may be made concerning the nonlinear 
orthotropic hardening coefficient. First, since the computations must 
be carried out at the Gauss point level, the computation of plastic 
material behavior will be expensive. The value of the coefficient is 
dependent on the anisotropy parameters, stresses and effective stress 
gradients at the end of the previous step, and the normalized stress 
Increment for the previous step. Note that it Is not a function of the 
magnitude of the prior step stress increment, but of the relative values 
of the stress increment components. In defense of the computational 

n p 

efficiency, the effective stress gradients rr 2 - are needed for later 

dS j 

use In the flow rule. The method is valid for orthotropic as well as 

Isotropic materials. For isotropic media, the method reduces to a 

von Mises criterion with hardening coefficients piecewise extrapolated 
-P 

along the s - e curve. The mechanics of implementation of the algorithm 
are presented in Appendix B. Appendix C shows the reduction of Hill's 
criterion to Mises' criterion for isotropic materials. 


Chapter 4 

APPLICATION TO COMPOSITES 

Application of the finite element technique discussed In Appendix 
A to laminated composites was previously discussed by Lin [25], 

Dana [26], and Dana and Barker [27]. Some of their conclusions will 
be repeated here for completeness. As discussed In Appendix A, the 
approach of the present study follows what has been termed the 
"macroscopic" approach. The approximations and Implications of this 
technique are discussed In the following paragraphs. 

4.1 Microscopic Versus Macroscopic Approaches 

The microscopic approach to the analysis of composite materials 
considers the fiber and matrix materials to be separate, each having 
the properties of the bulk fiber and matrix materials, respectively. 
Material properties are experimentally determined or assumed for the 
constituents, which may be isotropic or anisotropic. Some degree of 
bonding between the fiber and matrix is assumed. The degree of bonding 
is chosen based on the subject of the investigation. The microscopic 
approach appears to be useful for the study of fiber/matrix Inter- 
actions and failure mechanisms. However, It Is apparent that the 
computer cost for micromechanical analysis of laminates with even a 
small number of plies would be prohibitive. Indeed, Foye [35] used 
a large number of finite elements to perform a two-dimensional 
analysis of one fiber and its surrounding matrix, even after a 
significant degree of symmetry was assumed. 
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The macroscopic approach used In the present study assumes that 
the fiber reinforced material Is, on the Individual ply level, a 
continuous, homogeneous, orthotropic material. Properties used are 
generally determined by tests on unidirectional materials. A descrip- 
tion of the tests necessary to determine the elastic (and Inelastic) 
properties of an orthotropic material Is given In the report by 
Brinson and Yeow [62]. When using the macroscopic approach, the 
analyst should realize that all fiber/matrix Interactions actually 
occur at the microscopic (actually molecular) level, and have not 
been considered. Possible failure mechanisms such as fiber/matrix 
shearing, matrix cracking, and fiber fracture cannot readily be 
studied. However, Important effects such as Interlaminar stresses, 
stress concentrations near holes, free edges, and corners are well 
treated. 

In summary, one must conclude that the microscopic approach Is a 
more accurate model of what really occurs in a composite material, 
but has computational complexity which renders It Impractical for all 
but those problems involving single fibers and the surrounding matrix. 
The macroscopic approach yields considerable valuable Information, 
and Is computationally efficient for many laminates. Problem size 
and complexity are limited by core storage and cost. The indicated 
approach Is to use the microscopic model for matrix/fiber interaction 
studies and the macroscopic model for problems Involving stress 
concentrations, end and edge effects, and interlaminar stresses. 

As pointed out by various investigators [6,8,9], end and edge effects 
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are manifested In a region known as the "boundary layer" which extends 
Inward from the free edge of a laminate to a distance approximately 
equal to the thickness of the laminate. Away from the boundary layer, 
the computationally simple lamination theory solution Is recovered 
by finite difference and finite element techniques. 

The present study uses the macroscopic approach and Its primary 
topic Is the determination of Interlaminar stresses near holes and 
free edges. Several material models have been considered. Linear 
and nonlinear thermal expansions, as well as temperature dependent 
plasticity, have been treated. 

4.2 Modeling Considerations 

Care must be taken when using the finite element method that the 
assumptions of the method are not violated. In addition, the analyst 
must insure that the finite element model (or grid) used Is an 
accurate model of the real structure. Boundary conditions must 
match, and the finite element model must be capable of undergoing 
any deformation which Is possible 1r, the real structure. The former 
requirement Involves the use of sufficient nodal points on the 
boundary for adequate modeling. The latter condition Involves 
placing sufficient elements In the region away from the boundary 
layer. When the number of elements Is increased, care must be taken 
that element aspect ratios are in the proper range. Poor aspect 
ratios can lead to erroneous results. In addition, the analyst 
must be sure that finite elements which are large In one direction 
do not join finite elements which are small in the same direction. 
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Such a condition, Involving flexible (large) element connected to 
stiff (small) elements, may result In the assembled stiffness matrix 
being Ill-conditioned or even quasi-singular. 

The properties of the element used In the present study have been 
fairly well documented by previous Investigators [25,26,27]. Lin [25] 
determined that the best shape for the element for modeling Isotropic 
materials was a rectangular parallelplped, square (a*b) In the 12 node 
plane, (see Fig. 4.1) and with an aspect ratio of three (3c*a) In the 
eight node plane. Acceptable stiffness matrices were generated for 
width to thickness (a/c) ratios up to 12. Dana and Barker [27] 
reported that three elements were needed through the thickness of each 
ply of a laminated composite fo“ adequate determination of the through- 
the-thlckness stress gradients. This finding is consistent with that 
of Wang and Crossman [16], who used a large number of elements through 
the thickness In order to see the gradients. Consideration of the 
boundary layer dictates the size of the elements near boundaries, since 
large gradients exist in the boundary layer. The need for three 
elements through each ply plus acceptable aspect ratios provides 
another guide for element size. 

The requirements noted above are for very accurate elements. In 
regions where gradients are low (i.e. the lamination theory region), 
or where accuracy of stress is not critical (i.e. the region away 
from a hole being studied), elements with less than optimum size and 
aspect ratio may be used. This is common practice in finite element 
analysis. Although <\i1s practice is questionable from a theoretical 
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standpoint, it Is often the only way a problem can be run due to 
computer cost or storage limitations. As long as the analyst realizes 
the regions of accuracy and inaccuracy, the technique is viable. 

4.3 Grid numbering and bandwidth considerations 

As mentioned in Appendix B, the bandwidth of the global stiffness 
matrix is very important. But it is only the mean and not the maximum 
bandwidth that is crucial when using the skyline column storage scheme 
[63]. In this scheme, the upper triangular portion of the global 
stiffness matrix is stored, column by column, as a vector. Only the 
terms up to and including the upper nonzero entry in a column are 
stored. Larger mean bandwidth thus has the double effect of 
increasing the core required to store the matrix as well as the 
number of operations required for triangularization and back substitu- 
tion. For example, consider the grid and numbering of Figures 4.2 
and 4.3. The grids shown are similar to those used for simulation of 
the off-axis tensile specimens reported in the next chapter. The 
loading is uniform extension in the X-direction. The boundary 
conditions are that the ends of the specimen retain the same size as 
the deformation proceeds. The number of degrees of freedom of the grid 
is 168. Thus, the full upper triangular portion of the matrix would 
require 14196 core locations for its storage. The grid, numbered as 
shown in Figure 4.2, has 168 degrees of freedom, a maximum half 
bandwidth of 150, and required 12900 core locations for its storage 
(not counting the diagonal pointer array). The same grid, numbered 




FIGURE 4.3 GRID AFTER BANDWIDTH MINIMIZATION 
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as shown In Figure 4.3, also has 168 degrees of freedom, but has a 
maximum half bandwidth of 72, and requires only 7284 storage locations. 
It Is apparent that minimization of bandwidth should always be done 
before any attempt Is made to run problems. 

The bandwidth of all finite element grids used In this study was 
minimized using a modification of the computer program BANDIT [64] 
and the Glbbs-Poole-Stockmeyer [65] minimization algorithm. Other 
algorithms [66,67] are also available for bandwidth minimization. The 
experience of the author In using BANDIT Indicates that the Glbbs- 
Poole-Stockmeyer method Is generally cheaper and yields better results 
than the Cuthill -McKee [66] method. 


Chapter 5 

RESULTS AND DISCUSSION 

This chapter details the results of a series of problems, dealing 
primarily with laminated composites, which were solved using the 
techniques derived in the previous chapters. Comparisons are made to 
experimental data as well as the predictions of other models. Some 
of the problems were solved using the CDC 6700 computer system at 
the U. S. Naval Surface Weapons Center, Dahlgren, Virginia, and the 
remainder were solved on the IBM S370/3033 at VPI&SU. The computer 
program NALCOM (Nonlinear Analysis of Laminated COMposites) developed 
for this study is described in Appendix B. A listing of the program 
is available from the author on request. 

5.1 Isotropic Material 

As a check of the elastic-plastic capabilities of the program, an 
isotropic material was analyzed. The material was a chromium-nickel 
steel for which the stress-strain curve is presented as Figure 11 of 
Reference 60. The Ramberg-Osgood coefficients are given in Table 5.1. 
The yield stresses were chosen as the points of significant deviation 
from linearity. A comparison of experimental and computed values is 
shown in Figure 5.1. The program predicts linear elastic unloading 
and subsequent yield when the stress state reattains the yield 
hypersurface through some loading cycle. Also note that the isotropic 
hardening model does not predict any Bauschinger effect. 

As stated in the initial derivation, this is a known shortcoming of 
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o x VS e x FOR ISOTROPIC MATERIAL 


FIGURE 5.1 
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the Isotropic hardening model. 

5.2 Off-axis Tensile Specimens 

One common test for orthotropic materials which has received 
considerable attention Is the off-axis tension test [68,69 70]. Shown 
schematically In Figure 5.2, the test Involves applying tensile 
loading to a unidirectional material which Is oriented In the test 
fixture such that the fibers are at some angle e with respect to the 
load axis. Chamls [68] has proposed the use of a 10 degree off-axis 
tensile test for determination of the shear modulus 6^* Cole and 
Pipes [69,71] performed a series of off-axis tensile tests on NARMCO 
5505 Boron/Epoxy. Sandhu [70] found their 15° off-axis data suspect 
and performed the test again, with significantly different results. 

The computer program NALCOM was used to predict the linear elastic 
and elastic-plastic response of Boron/Epoxy off-axis tensile specimens. 
Material properties are shown in Table 5.1. The finite element grid 
used for all off-axis tensile modeling is shown in Figure 5.3. The 
model is symmetric about the x-y plane. The ends are assumed to remain 
straight and are free to move only in the x-direction for the linear 
elastic computations. For the nonlinear predictions, the load was a 
uniform stress applied to the end of maximum x, which was constrained 
to no rotation and no y or z direction motion. These end conditions 
were intended to simulate perfect gripping, no end tabs, and no 
rotation of the grips. 

5.2.1 Linear Elastic Off-axis Tensile Specimen Modeling 

Off-axis tensile specimens of NARMCO 5505 Boron/Epoxy and 
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T300/5208 Graphite/ Epoxy were modeled using linear elastic material 
properties from Cole and Pipes [69,71] and unpublished data obtained 
from the Lockheed Company. Unless otherwise specified, all experi- 
mental data for T300/5208 were taken from the Lockheed data package. 
The finite element grid shown In Figure 5.3 was used for these 
analyses. The dimensions are given In Table 5.2. To simulate the 
effect of bonding end tabs to the composite as an aid to gripping, 
a length at each end equal to 20% of the total coupon length was 
constrained to no y-dlrectlon motion. This has the same effect as 
a perfectly bonded tab which allows extenslonal but no In-plane 
bending deflection along the tab length. 

The linear elastic characteristics of Boron/Epoxy off-axis 
tensile coupons as determined by an x-dlrectlon strain loading are 
shown In Figures 5.4 and 5.5 and Table 5.2. The quantity S^j as 
computed from the transformation equations Is Included for comparison. 
Note the close agreement between the model and Cole and Pipes' 
experimental data with the exception of the apparent modulus E xx at 
15°. As previously mentioned, Sandhu [70] found the 15° data 
suspect. Sandhu 's 15° data is also shown on Figure 5.4, with 
significantly improved comparison. Note the small effect of the 
end tabs (Table 5.2). The effect of the end tabs is to decrease 
the coupon aspect ratio L/h. Since Cole and Pipes used aspect ratios 
ranging from 18.7 to 28, even the decrease due to the tabs is not 
sufficient in most cases to bring the ratios below the value of 
14 to 16 which Is generally regarded as an acceptable value for 
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FIBER ANGLE u (DEGREES) 

FIGURE 5.4 LOAD DIRECTION MODULUS VS, FIBER ANGLE 
FOR OFF -AX J S TENSION OF BORON/EPOXY 





FIGURE 5.5 MAJOR POISSON'S RATIO VS. FIBER ANGLE 
FOR OFF-AXIS TENSION OF BORON/EPOXY 
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elimination of end effects. Figure 5.6 is a plot of the deformed 
finite element grid for Boron/Epoxy and o * 15°. Note the In-plane 
warping which Is characteristic of off-axis tensile specimens. Such 
warping gives rise to substantial variations In stress and strain 
across the width of the specimen. Such variations raise the question 
of accuracy of strain measurements on such specimens, since strain 
gauges are of finite dimensions. 

The computed linear elastic moduli and Poisson's ratios of 
Graphite/ Epoxy are shown in Figures 5.7 and 5.8. No experimental data 
appear to exist for off-axis tensile testing of Graphite/Epoxy, so no 
comparisons are possible. The basic shape of the curves is the 
same as for Boron/Epoxy. Note that the major Poisson's ratio is 
again a maximum at approximately 20°. 

5.2.2 Nonlinear Response of Off-axis Tensile Specimens 

As a check of the orthotropic plasticity model developed in 
Chapter 3, the elastic-plastic response of Boron/Epoxy off-axis 
tensile specimens was computed for comparison to the nonlinear data 
of Cole and Pipes [69]. As with the linear case, the 15° data of 
Sandhu [70] was also considered. Ramberg-Osgood coefficients were 
determined from the 0°, 90°, and pure shear (from tube torsion) data 
of Cole and Pipes. Comparisons for Boron/Epoxy were made to two 
sets of experimental data and two other analytical models. The 
elastic-plastic computations were also run for Gr/Ep and are 
presented without comparison, since no off-axis tensile data 
for Gr/Ep are known to exist. Comparisons of all the models to 



FIGURE 5.6 DEFORMED GRID OF 15° OFF-AXIS TENSILE SPECIMEN 
(f/h = 23.3, e x = 0.1) 
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FIBER ANGLE 6 (DEGREES) 


FIGURE 5.7 LOAD DIRECTION MODULUS VS. FIBER ANGLE 
FOR OFF-AXIS TENSION OF GRAPHITE/EPOXY 
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FIGURE 5.8 MAJOR POISSON'S RATIO VS. FIBER ANGLE 
FOR OFF-AXIS TENSION OF GRAPHITE/EPOXY 
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the Boron/Epoxy data are good, with the results of the present study 
showing better agreement In some areas and worse In others. Note that 
the present study deviates from experiment primarily when plastic 
she'.r strains are more significant. Plastic shear strains were 
determined to exist at all angles other than 0° and 90°. However, 
at 15°, 30°, and 45°, the plastic shear straining was noted to be the 
dominant plastic strain component. Since agreement between computed 
and experimental results is good for angles where plastic shearing 
strain is relatively low and significant differences occur when plastic 
shear strains are large, one must question either the manner in which 
plastic shear strains are included or the numerical values which were 
used as input. It appears that for all angles the plastic strain is 
initially underestimated, gradually becoming overestimated. The 
program NALCOM allows a two segement Ramberg-Osgood approximation 
to the stress-strain data. The method of determination of the 
coefficients is to plot stress versus plastic strain on a log-log 
scale and determine the best bilinear or linear fit. The tube 
torsion data from Tube 11 1 B from Cole and Pipes [69] is plotted in 
this manner in Figure 5.9. It appears that the data are not well 
approximated by a bilinear fit. This assures that the model will 
have difficulty in duplicating the experimental data. Fortunately, 
a linear or bilinear Ramberg-Osgood fit is sufficient for most 
engineering materials. 

The nonlinear response of Boron/Epoxy off-axis tensile specimens 
is shown In Figures 5.10-5.15. The results of the present study are 
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compared to the experimental results of Cole and Pipes [69] for all 
angles, the experimental results of Sandhu [70] for 15°, the computed 
results of the model of Hashin [12,13] as programmed by Plndera [72], 
and a modified form of Hashln's model due to Plndera [72]. 

The 15° off-axis response of Boron/Epoxy Is shown In Figure 5.10. 
The shearing effects are strong at 15°, and It Is believed that this 
fact in conjunction with the previously mentioned inability to 
accurately prescribe the shear behavior with a two segment Ramberg- 
Osgood approximation is at least part of the reason for the difference 
in computed and experimental results. Also, the Ramberg-Osgood 
parameters were taken from Cole and Pipes' data, which show consider- 
able variations from one specimen to the next. Recall that the 15° 
data of Cole and Pipes is suspect, adding additional complication to 
meaningful comparisons. 

The off-axis response of Boron/Epoxy at angles of 30°, 45°, 60°, 
75°, and 90° is shown in Figures 5.11, 5.12, 5.13, 5.14 and 5.15, 
respectively. The comparison at 90° is very good, and the discrepancy 
at the other angles is again believed due to the poor shear data used 
as input to NALCOM. 

The nonlinear off-axis response of Gr/Ep is shown in Figure 
5.16. As with the linear elastic results, no experimental results 
are available except at 0°, which is linear. The angles 15°-75°, 
in which plastic shear straining is large (i.e. 15°-45°) show signi- 
ficant nonlinearity and become increasingly linear as the linear 
behavior of the 90° direction is approached. Note that the response 
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FIGURE 5.11 OFF-AXIS TENSION OF BORON/EPOXY, FIBER ANGLE = 30° 
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FIGURE 5.12 OFF-AXIS TENSION OF BORON/EPOXY, FIBER ANGLE « 45° 
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FIGURE 5.13 OFF-AXIS TENSION OF BORON/EPOXY, FIBER ANGLE = 60° 
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FIGURE 5.14 OFF-AXIS TENSION OF BORON/EPOXY, FIBER ANGLE « 75° 
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of e = 90° is linear to failure. 


5.3 Laminate Studies 

The final test of the model developed for the present study is 
to compute the residual stress state resulting from the curing of 
symmetric composite laminates. Various geometries and material 
models were considered. Laminates were limited to symmetric layups 
consisting of two fiber orientations. This limitation was introduced 
in order to hold the expenditure of computer resources to a reasonable 
level while considering problems of significance. Comparisons to 
experimental data and the results of other models are made whenever 
possible. Such comparisons are often complicated for a variety of 
reasons. For example, Rybicki and Schmuesser [32], Pagano [6], and 
Wang and Crossman [16,73] used material models which set all three 
shear moduli (G-^. G 13’ G 23 J and all three Poisson's ratios 
(v 23 » v i 3 » v i2^ ec l ua ^ • These assumptions are inconsistent with the 
theory of elasticity for transversely isotropic materials and also 
with the material model of NALCOM. In addition, the models of Wang 
and Crossman [73] and Renieri and Herakovich [40] are quasi -three 
dimensional, while NALCOM is fully three dimensional. The nonlinear 
material model used by Herakovich and his associates [40,74,75] is 
nonlinear elastic, while the material model imp’emented in NALCOM is 
based on a plasticity theory. The only other study known to have 
considered second order thermal effects is that of Hahn and Pagano [15], 
which was used in conjunction with a lamination theory model and 
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other simplifying assumptions. Dana [26] and Dana and Barker [27] did 
not consider thermal effects. In the light of these differences in the 
models, comparisons are often possible only in a qualitative sense. 

The basic laminate geometries for the present study are shown in 
Figures 5.17 and 5.18. For simplicity, only one-eighth of the laminate 
is shown. For symmetric cross-ply (i.e. e =0,90) laminates, there is 
symnetry about all three coordinate planes, and one eighth of the 
plate is sufficient for computations, provided that the loads are 
also symmetric. For symmetric angle ply laminates, there is symmetry 
only about the x-z plane, so one-half of the laminate must be modeled. 
As previously stated, three elements are required through the thickness 
of each fiber orientation. The object of the present study is to 
examine the stress fields near the edge and end of the solid laminates 
and near the hole in the laminates with central holes. Accordingly, 
the finite element grids used for the solid and central hole cross-ply 
laminates are shown in Figures 5.19 and 5.20, respectively. For 
accuracy in angle ply laminates, the upper half of the laminates must 
be modeled. This requirement will be discussed later in more detail. 

The one eighth symmetric sections (Figures 5.19 and 5.20) contain 
six layers of four elements each. More elements in each layer would 
have yielded better results, but at considerable expense in computer 
time and storage. Each element contains 32 (4x4x2) Gauss points. 

The stress is computed at each Gauss point and averaged through the 
element thickness to obtain through the thickness stress distribu- 
tions. Of primary concern is o z> the interlaminar normal stress. 






FIGURE 5.18 GEOMETRY OF ONE-EIGHTH OF A LAMINATE WITH CENTRAL HOLE 










FIGURE 5.19 FINITE ELEMENT GRID FOR LAMINATE EDGE STUDIES 
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Results are computed for both linear elastic and elastic-plastic 
material models with constant, first order nonlinear, and second 
order nonlinear thermal behavior. The following paragraphs detail 
these analyses. 

5.3.1 Uniform [0/90] s Laminate Studies 

5. 3. 1.1 Strain Loading 

As a first comparison to other models, a [0/90] s boron/epoxy 
laminate was subjected to a uniaxial strain loading. Figures 5.21 
and 5.22 show the comparison of results of the present study with the 
results of Pagano [6] and Wang and Crossman [73]. Note that the 
comparison is excellent. The material model of NALCOM was modified 
to conform to the simplified model of the other investigators for 
this comparison. 

5. 3. 1.2 Thermal Loading 

T he next problems of concern were the linear elastic cooling from 
the elevated curing temperature to room temperature followed by 
imposition of mechanical load. For these analyses, the three thermal 
effect options (i.e constant properties, first order nonlinear and 
second order nonlinear) were considered. The finite element grid of 
Figure 5.19 was subjected to temperature changes totaling -180°F, 
in one step for the constant property case, using average properties 
over the range, and in nine steps of -20°F each for the nonlinear 
cases. The material for these analyses was T300/5208 graphite/epoxy. 

The material properties used are given in Table 5.3. As previously 
stated, the interlaminar stresses are of primary concern for this study. 
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FIGURE 5.21 TRANSVERSE DISPLACEMENT ACROSS TOP SURFACE (z = 2h), [0/90] B/Ep LAMINATE, 
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TABLE. 5.3. ELEVATED TEMPERATURE PROPERTIES OF T300/5208 GRAPHITE/ EPOXY 
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Figure 5.23 shows residual Interlaminar normal (o z ) stress vs. 
y/b near the midplane (z/2h * 0.035) of a [0/90] s laminate of 
Graphite/ Epoxy. Unless otherwise stated, all laminate analyses were 
performed for T300/5208 graphite/epoxy. Figure 5.23 shows the 
effect of the three thermal models. The effects are small, amounting 
to maximum differences of approximately 10%. Note that in the 
region away from the boundary layer (y/b < 0.6) the lamination theory 
solution (o z * 0) Is recovered exactly. An earlier version of 
NALC0M with no equilibrium iteration did not recover the lamination 
theory solution at y/b * 0. The curves of 5.23 were generated using 
the elastic-plastic analysis of NALC0M. The maximum value of Hill's 
yield function at z/2h - 0.035 was 0.419. This occurred In the region 
where the end boundary layer interacts with the edge boundary layer 
(x/a * 0.65, y/b = 0.89). At this point It should be emphasized that 
the three dimensional analysis shows the existence of the end (x/a « 1) 
boundary layer as well as the edge (y/b = 1) boundary layer normally 
discussed by Investigators using two-dimensional or quasi three- 
dimensional models. The boundary layers Interact at corners (x/a > 0.6, 
y/b > 0.6), often with surprising results, which are dependent on 
through the thickness location and the laminate layup. For linear 
elastic models, the interaction Is basically a superposition, but the 
Introduction of a yield condition and flow rule results in a 
considerably more complex situation. 

Figure 5.24 shows residual o z vs y/b at x/a = 0.093 at the 0/90° 
interface for the [0/90] s Gr/Ep laminate. The difference between 
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the different thermal models Is again approximately ten percent. The 
lamination theory solution of a z = 0 Is again satisfied for y/b < 0.6. 

As before, previous results with no equilibrium Iteration did not 
satisfy this condition. Thermal loading was accomplished incrementally, 
and equilibrium was normally achieved in one to two Iterations at each 
Increment. Yielding of the material on both sides of the 0/90 
Interface was predicted. The amount of plastic deformation was small, 
and stress differentials between elastic and elastic-plastic models 
is thus small. Yielding occurred near the corner (x/a * 1, y/b * 1) 
of tie laminate, with yielding being more pronounced in the top (0°) 
layer. The reason for yielding at the 0/90 interface when none 
occurred at the midplane is the thermal expansion mismatch of the 
0° and 90° directions, resulting in large axial ( a ) and transverse 
( 0 ^) stresses near the corner. Since yielding depends on absolute 
value of stress, even though the algebraic sign of o and a are 

X y 

different, the combination results in yielding. Values of interlaminar 
shear stress (t xz and T yz ) are also large near the corner. This is a 
good example of how a three-dimensional analysis reveals phenomena 
not treated by two dimensional models. No yielding was predicted for 
x/a = 0 (corresponding to a 2-D analysis of [0/90] $ ) or y/b = 0 
(corresponding to a 2-D analysis of [90/0] s ), but the interaction was 
sufficient to result in yielding. Yielding occurred in the -20°F 
thermal increment from 90°F to 70°F. 

Figure 5.25 shows the variation of o z through the laminate 
thickness near the edge (y/b = 0.977) and away from the end 



FIGURE 5.25 a 2 VS. z/2h, [0/90] s Gr/Ep LAMINATE (y/b = .977, x/a = .043) 
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(x/a = 0.043). The effect of the different thermal analyses is more 
pronounced near the midplane, with little difference in the upper 
(0°) layer. All models converge to the equilibrium solution (o z = 0) 
at the top surface (z/2h * 1 . ) . 

Several comments concerning the laminate analyses may be made. 

The addition of second order nonlinear thermal effects was found to 
have the effect of reducing the axial (a x ) stress near the origin 
(x/a = y/b = 0), as predicted by Hahn and Pagano [15]. This effect 
cannot be ignored, especially in the investigation of elastic-plastic 
interlaminar effects. When yielding is imminent, even small changes 
in stress become significant. 

5. 3. 1.3 Combined Mechanical and Thermal Loading 

All thermal cooling (aT = -180°F) was followed by the application 
of a uniform x-direction traction (o^) of 1000 psi. Although the 
magnitude of the load was small, the effects are significant. All 
the Gauss points near the free corner which were exhibiting plastic 
behavior started to unload in the plastic sense. This is easily 
explained if we consider that certain of the stress terms exhibit 
reversal in algebraic sign when going from thermal cooling to 
x-direction extension. Had the x-direction traction been continued, 
yielding would undoubtedly have been resumed, although not necessarily 
in the same location. Again, the three-dimensional effect of the free 
corner must be considered. The two dimensional and quasi three 
dimensional analyses of previous investigators is valid away from 
three dimensional effects of ends and holes, but one must realize the 
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limitations. Also, one must be aware that Intuition gained from examin- 
ation of two dimensional analyses may not extend to three dimensional 
problems, and may actually be misleading If used In attempts to predict 
three dimensional effects. 

5.3.2 Uniform [+e] $ Laminates 

Also of concern Is the response of symmetric angle ply laminates. 
Previous investigators [26,32,40] have analyzed [±e] $ configurations 
( 6^0,90) assuming varying degrees of synmetry. With the exception of 
Dana, et al [26,27], no consideration seems to have been given to the 
question of symmetry. In order for symmetry to exist about a plane, 
all nodal points originally on the plane must undergo displacements only 
in that plane . A synmetry plane may be thought of as a mirror, where 
all quantities (material properties, displacements, loads, stresses) 
must be mirrored by the plane. If we consider a [±e] $ laminate, 
e^0°,90°, the only symmetry which exists is across the midplane (Z=0). 

The symmetry conditions for a long laminate advanced by Hsu and 
Herakovich [8] are 


u(x,y,z) 

= u(x,y,-z) 

(5.1a) 

v(x,y,z) 

= v(x,y,-z) 

(5.1b) 

w(x,y,z) 

= -w(x,y,-z) 

(5.1c) 

v(x,y,z) 

= -v(x,-y,z) 

(5. Id) 

w(x,y,z) 

' w(x,-y,z) 

( 5 . 1 e ) 


plus an experimentally verified condition [8] that 


u(0,y,h) = -u(0,-y,h) 


(5.2) 
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where Z=h Is the top surface. The condition (5.2) Is generalized to 

u(0, y,z) = -u(0,-y,z). (5.3) 

NALCOM was used to model the top half of a [±45] s laminate, with one 
element being used for each quadrant of the laminate, and six elements 
through the thickness as before. The load was a uniform x-dlrectlon 
extension of e =0.001. For uniform x-dlrectlon extension, if symmetry 

A 

exists about the x-z plane as proposed by other Investigators, the 
y-directlon displacements along the Initial lines of y=0 should all 
be zero. Flgurpr 5.26 and 5.27 show the seven lines of Initial y=0 
after load application for thick (a/h=15) end thin (a/h=150) laminates, 
respectively. Note that the only locations where all are zero Is at 
the ends (x/a=±l) where the ends were rigidly clamped, and at the 
axial center, x/a=0. Everywhere else the displacements are small, 
but nonzero. The symmetry condition is thus seen to be nonexistent. 

The assumption of symmetry about y=0 is seen to be increasingly poor 
as the laminate becomes thicker. The coarse grid used for these 
analyses did not permit what were considered accurate stress 
determinations. It is possible that while the violation of the 
y=0 symmetry may be small with regard to displacements, the symmetry 
of stresses which results may contain severe inaccuracies. The 
symmetry conditions which are not satisfied are 5. Id and 5 . 1 e . For 
infinitely long laminates, the conclusion is that the usual symmetry 
assumptions are valid. However, three dimensional studies (e.g. [30]) 
which have assumed symmetry must be questioned, especially in regions 


, z/2h)/2b X 10' 
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near planes where symmetry has been Improperly assumed. By Invoking 
St. Venant's Principle, one might be able to use results from other 
regions of the laminate. Additional Investigation Into this problem 
Is Indicated. 

5.3.3 Laminates with Central Circular Holes 

Due to core limitations, only [0/90] s laminates with central 
circular holes were modeled. The geometry Is shown schematically In 
Figure 5.18, and the finite element grid Is shown In Figure 5.20. 

The material system considered was T300/5208 graphite/epoxy. As for 
the uniform laminates, the Interlaminar normal stress (o z ) was of 
prime concern. 

Figure 5.28 shows c> 2 vs. position around the hole near the mid- 
plane. The difference In the three thermal models Is more pronounced 
at the extreme values of «j>. The difference Is approximately ten 
percent. No yielding of the material was predicted at the midplane. 
Figure 5.29 shows o z around the hole at the 0°/90° interface. 

Yielding was initiated In the upper (0°) layer near 6*45°. The degree 
of plasticity was small, resulting In little stress relief. Yielding 
occurred during the -20°F Increment from 90°F to 70°F. As with the 
uniform laminates, the thermal cooling was followed by a uniform x- 
dlrectlon traction of 1000 psl. The yielded a**ea unloaded elastically, 
the value of Hill's yield function undergoing a slight decrease. 
Interlaminar shear stresses t and t were significant ("1400 psl) 

Jr fc At 

near the Interface and e=45°. 


Figures 5.30 and 5.31 show the through the thickness variation 
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of o 2 near the hole for ♦■S.U and 86.87°, respectively. The effect 
of the different thermal models Is again approximately ten percent, 
with the differences more pronounced near the midplane. The 
equilibrium condition c? 2 «0 Is recovered for z/2h«l. 


Chapter 6 

SUMMARY AND CONCLUSIONS 

The present study was concerned with a number of facets of the 
analysis of laminated composite materials. A plasticity theory, based 
on Hill's orthotropic yield criterion and the Incremental flow theory 
of plasticity was formulated and made compatible with the finite 
element method of stress analysis. Techniques of Including nonlinear 
hardening of orthotropic materials via unidirectional Ramberg-Osgood 
coefficients as well as temperature dependent plasticity, and linear, 
firsc and second order nonlinear thermal expansions was included. 

A computer implementation was formulated, coded, and utilized to study 
a number of problems relevant to composite materials. 

Based on the work done in this study, the following conclusions, 
comments, and recommendations for future research are presented. 

1. Three dimensional effects are definitely present In 
laminate composites, especially near free edges, holes, and corners. 
Treatment of these effects requires a fully three dimensional analysis 
capability. 

2. For problems (e.g. cross-ply laminates, long laminates) 
which are closely approximated by generalized plane strain models, 
quasi-three dimensional analyses agree well with fully three dimension- 
al results. They have the advantage of substantially lower cost, 

plus easier Input generation and output assimilation. 

3. Nonlinear behavior is a factor in compos ite materials, and 
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can significantly alter stress distributions. However, plasticity 
theories for orthotropic materials are not well developed, and the 
analyst must take care that spurious results do not lead to wrong 
decisions and faulty structural designs. 

4. In material systems which are temperature dependent, as are 
composites, nonlinear thermal effects may be Important, and should 
not be Ignored. Where stress and/or temperature gradients are high, 
second order nonlinear thermal effects may also be Important. 

5. Symmetry conditions which have been advanced and widely 
used for angle ply laminates appear to be Identically correct for an 
extremely limited number of positions In the laminate. Both the 
analyst and experlmentor should take care not to assume symmetry where 
none exists. In addition, experimenters should be continually aware 
of the large gradients of strain which may exist when composites are 
extended off-axis. Variations may exist over the length of a strain 
gauge. 

6. The initial strain technique of finite element analysis Is 
a viable technique for modeling thermal and plasticity effects. It 
has the severe drawback of not being able to treat properties which 
differ in tension and compression. It Is possible that a tangent 
stiffness method could be made competitive with respect to execution 
time and have the advantage of different properties In tension and 
compression. 

7. An element with as many degree of freedom as the one used 
in the present study rapidly lead to stiffness matrices with large 
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bandwidth, often resulting In catastrophic storage problems. Out of 
core solvers are very slow In assembly and suffer somewhat In 
caution. A minimization technique which does not require assembly 
of the global stiffness matrix is Indicated as a viable solution 
algorithm. 

8. Thick cross-ply graphite/epoxy laminates are not highly 
susceptible to plasticity during the cooling cycle, even when free 
edge influence is considered. Were the transverse modulus raised, 
plasticity might be a problem. 

9. When modeling nonlinear effects, an equilibrium iteration is 
essential. Omission of equilibrium iteration or failure to converge 
can result in instabilities which render useless any values computed 
after deviation from equilibrium. 

10. When using any nonlinear model and describing stress-strain 
relationships, the tangent modulus is of the utmost importance. Even 
if a is correctly computed given e, a faulty tangent modulus yields 
poor results when used to extrapolate to the next step. Exponential 
models such as Ramberg-Osgood suffer severely from this problem since 
small deviations raised to a power greater than unity are magnified. 
Such errors propagate and grow very rapidly. 
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APPENDIX A 

LINEAR ELASTIC FINITE ELEMENT RELATIONSHIPS 

The finite element used Is shown In Figure A.l. The sides with 
four nodes are described by the cubic "serendipity" Interpolation 
functions as described by Zlenklewlcz [76]. The details of the Isopara 
metric finite element concept are also given In Reference 76. The 
technique basically Involves mapping a distorted shape In the 
Cartesian (x,y,z) coordinate system Into a cube In a local U, n,0 
coordinate system where L, n, and c range from -1 to +1. The 
relationship between the global Cartesian and the local curvilinear 
coordinates Is 

x * + NgXg + ... + N 24 x 24 * N^ 1 * 1,24 

y * N 1 y 1 + N^ + ... + N 24 y 24 * N^ 1-1 .24 (A.l) 

z - N 1 z ] + + ... + N 24 z 24 • N 1 z i 1 » 1,24 

where the N^ are the Interpolation functions for the 24 nodal points 
and the x^, y^, z^ are the Cartesian coordinates of the nodes. If 
we Introduce the notation 


8 ^i n 0 “ nn 1 c 0 “ w 1* 
then for the corner nodes c * ±1 * ±1 ^ ■ ±1 


(A. 2) 


N 1 8 87< 1+ 5 0 )( 1+ n 0 )(l + t 0 )[9(c 2 +n 2 )-10] 


1 « 1 ,4,5,8,17,20,21 ,24. 


(A. 3) 



Ill 


For the triside nodes * ±1 Hj • t | “ tl 

N 1 " CT< 1+ «o )(U9n O )(1n O )(1 ‘ n2) 1 " 9,10.11 .12,13,14,15,16. 

(A.4) 

For the triside nodes * *1 ^ ■ t y ■ ±1 


N 1 * |y< 1+9 e 0 )(1+n 0 )(ln 0 )(1-C 2 ) 1 » 2,3,6,7.18,19,22,23. 


(A. 5) 


For an Isoparametric element, the same Interpolation functions are 
used for the assumed dlsp’acements as for the geometry. Hence, 


u ■ 

v ■ N,*, 
w ■ N^w,j 


(A. 6) 


where u,v, and w are the x,y, and z displacements, respectively, and 
the u^, v^, and w^ are the values at the 1^ node. 

The strain-displacement relationships are derived based on small 
strain-small displacement theory. For the three dimensional case, 
these relationships may be written as 


yz 


'xz 


xy 


3 

dx 


!=■ o 


3_ 

3y 


3_ 

3Z 


3_ 

9Z 

3_ 

3y 


3 

3X 


3_ 

3Z 

3_ 

9y 


0 h 


0 


(A. 7) 
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Substituting (A.6) Into (A. 7), 



where 

I Is the 3x3 Identity matrix 

(q) Is the 72x1 vector of nodal displacements given by 



We define the [B] matrix (strain-displacement relationships) by 


{c } - 


(A. 10) 
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By comparison of (A. 10) and (A. 8), 


[B] - 


3_ 

3X 

0 

0 

0 

3 

sy 

0 

0 

0 

3_ 

3Z 

0 

3_ 

3Z 

3 

sy 

3 

3Z 

0 

3_ 

3X 

3 

sy 

3 

3X 

0 


[IN 1 |IN 2 |,..|IN 24 ]. (A.ll) 


From (A. 3), (A. 4), and (A. 5), we recall that the are functions of the 
local coordinates 5 , n» 5. In order to determine the elements of [B] 
we require a relationship between the derivatives in the global 
(x,y,z) and local U.n.t) coordinate systems. This relationship is 
given by the Jacobian matrix [J] of the transformation where 


[J] = 


3X 3y. 3 z 

35 35 35 

3X_ 3y. 3Z_ 

3n 3n 3n 

3x 3y 3z 

35 35 35 


(A. 12) 


= [J] 


(A.12a) 
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Inverting (A. 12a), 



(A. 13) 


where [J]~^ is the inverse of [J]. Note that if [0] is singular, the 
inverse does not exist. Physically, this indicates that the 
x,y«z-*£,n,c mapping is improper. Furthermore, if |J| (determinant of 
[J]) is negative an improper mapping is indicated. Improper mappings 
indicate faulty geometry input or elements which are so distorted 
that a unique mapping is Impossible. Substituting (A.l) into (A.12), 


[J] 


3N ] 

aT 

9N 1 

an - 

3N ] 


9N 


24 


35 

f?24 

an 

3N 


24 


at 



x 24 y 24 z 24 


(A. 14) 


For the application of finite elements to laminates, several 
assumptions are made. The material properties are assumed to be 
transversely isotropic, with the stiffer of the two directions always 
being in the plane of the element. The fiber direction so defined is 
oriented at an arbitrary angle. The material is assumed to be 
homogeneous. This approach has been termed the macroscopic approach, 
as opposed to the micromechanics approach of considering the fiber 
and matrix materials individually. In reality, one cannot ignore the 
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fiber-matrix Interactions which exist on the microscopic level. In 
practice, microscopic analyses cannot be handled for problems of 
significant size. The purpose of this Investigation Is to determine 
the stress distributions between laminae, an effect which Is 
manifested near the microscopic level. However, the physical 
dimensions of the solids to be treated preclude a microscopic analysis. 
Macroscopic analysis techniques will therefore be used with the tacit 
understanding that microscopic effects of possibly significant 
magnitude have been omitted. 

The elasticity matrix [D] for a homogeneous orthotropic 
material Is r -i 


[D] = 


where 

D 11 


'11 


1 ” v 23 v 32 


D .1^31 E 
u 22 F fc 2 


°44 " 6 23 


°12 °13 0 

0 

°22 °23 0 

0 

°33 0 

0 

°44 

0 

symmetric 

°55 

v 12 tv 13 v 3? , 

12 ' F l 22 

°13 

v 23 +v 21 v 13 r 

n 

23 ~ F t 33 

°33 

55 = G 13 

°66 


0 

0 

0 

0 

0 

3 66 


(A. 15) 


v 13 +v 12 v 23 


■33 


] ” v 12 v 21 


■33 


F = ^ " v ! 2 V 21 ~ v ! 3 V 31 ” v 23 v 32 ~ V 1 2 V 23 V 31 ’ v 21 V 1 3 V 32 ' 


The stress-strain relations for the material in the 
axis system (1,2,3) Is 


principal matrrlal 
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(A. 16) 


where the 1-direction is assumed to be coincident with the fibers. 
The material is isotropic in the 2-3 plane, leading to the reduction 
of the nine material constants of (A. 15) to the five required for a 
transversely isotropic material. The additional relationships due 
to the 2-3 isotropy are 


E 


2 



6 31 * 6 12 V 21 = V 31 

E 2 

G 23 = 2(1+v 23 ) • 


(A. 17) 


If the in-plane principal material axes (1,2) are not coincident with 
the in-plane global axes (x,y), the equations (A. 16) must be modified 
to result in a stress-strain relationship in the global system. Since 
and are actually second order Cartesian tensors and follow the 
transformation law of second order tensors. 


( 1 ) 

kl 

a a e (1) 
a ik a jl e kl ’ 


°1J = a ik a jl°kl 


'1j 


(A. 18) 


where and e^j are referred to the global axes, and are 


i j 


referred to the material axes, and a.^ are the direction cosines of the 
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global axes with respect to the material axes. For an angle e 
measured from the x axis to the 1 axis. 



cose -sine 0 

sine cose 0 , 

0 0 1 


(A. 19) 


which corresponds to an angular rotation e about the 3 (coincident 
with z) axis. Performing the operations (A. 18) with (A. 19), and 
replacing tensor strain e.. (i^j) with engineering strain y. . (i^j), 

I J I J 



(A. 20) 


where 

[t 2 ] 


m 2 n 2 0 

n 2 m 2 0 

0 0 1 

0 0 0 

0 0 0 

-2mn 2mn 0 


0 0 mn 

0 0 -mn 

0 0 0 

m -n 0 

n m 0 

0 0 m 2 -n 2 


(A. 21 ) 


and m = cose n = sine. 

Furthermore, the stresses in the global and material axis systems are 
related by 
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(A. 22) 


where 


[T,] = 


m 2 n 2 0 0 

n 2 m 2 0 0 


0 0 10 

0 0 0 m 

0 0 0 n 

-mn mn 0 0 


0 -2mn 


n m 2 „2 
0 m -n 


(A. 23) 


and m = cose, n = sine. The matrices [T^] and [Tg] are referred to 
as the stress and strain transformation matrices, respectively. Note 
that the presence of the factor two in [T^ and [T 2 ] In the shear 
terms accounts for the difference in tensor shear strain and engineering 
shear strain. 

Substituting (A. 21) and (A. 22) into (A. 16) and premultiplying both 
sides by [T ]""* , 
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^ E x | 



| e y I 


} - cv’mtv < 

|V | 

M 

X 

fr - 


f Y xz 1 



Y xy j 


Equations (A. 24) are the stress strain relationships in the global 
coordinate system. The rotated elasticity matrix [D] is defined by 

[D] = CT 1 ]* 1 CD][T 2 ]. (A. 25) 

The finite element used is formulated based on linear elasticity 
and the theory of minimum potential energy. Material nonlinearities 
may be introduced as pseudo-loads acting on the elastic body. This 
technique is discussed in Chapter 3. 

The total potential energy, n, of a given finite element is the 
sum of the strain energy, U, and the work of the external loads, V. 

In the following developments the matrix notation will be dropped and 
index notation used. The strain energy 1) of the element is 



si 


j E ij dn 


(A. 26) 


where 

is the stress tensor in globai (x,y,z) axes 
is the elastic strain tensor in global (x,y,z) axes 
n is the volume of the element. 

If we Introduce the initial strain tensor c?., (A. 26) becomes 
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u ■ 7 ! a (A ' 27) 

where is the total strain tensor. The initial strains may be 
due to thermal expansion, 

e tj = v (A - 28) 

where 

is the matrix of thermal expansion coefficients 
e is the temperature change from some reference state. 

Initial strains may also be due to moisture expansion 

•ij - v < A - 29 > 

where 

is the matrix of moisture expansion coefficients 
y is the moisture concentration. 

Initial strains may also be used to include material nonllnearitles 
(either static or time dependent). The initial strain method for 
nonlinear analysis is discussed in the text by Desai and Abel [50]. 

With the inclusion of initial stresses, o^, the linear elastic 
constitutive equations become 


°1j 


Vij * °ij 


(A. 30) 


where D^j is as defined by (A. 25) and the initial stresses may be a 
residual stress field due to some previous loading. Inserting (A. 30) 
into (A. 27), 


u = 7 f n C(e kj' E kj )D 1k (e 1j' e V + 2o ij (e 1j’ e ij ):Jdn 


(A. 31) 
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and 

u " \ f n tc M D 1k € 1j +e kj D 1k c 1j‘ 26 kj D 1k c 1j +2o 1j E 1j" 2o 1j e 1j ]dfi ^ A,32) 
The potential of the external loads, V, Is 

V « -q i Q i (A. 33) 

where 

Is the vector of nodal displacements 
Is the applied mechanical load vector. 

The total potential energy, n, of the element Is 

n = U + V (A. 34) 

where U and V are given by (A. 32) and (A. 33), respectively. Converting 
(A. 32) and (A. 33) to vector-matrix notation and inserting into (A. 34), 

n * } / n [L E j[D]{e}+Le 0 J[D]{ e °}-2LEj[D]{e 0 }+2LeJ{a o :-2La 0 J{e 0 }]dn-LqJ{Q}. 

(A. 35) 

Using (A. 10) in (A. 35), 

n = \ f a [LqJ[B] T [D][B]{q}+Le°J[D]{E°}-2LqJ[B] T [D]{e 0 } 

+ 2LqJ[B] T {o°}-2Lo 0 J{c 0 }]dfi-LqJ{Q}. (A. 36) 

Taking the variation of n with respect to the unknown nodal displace- 
ments (q } and setting to zero, 

6n * 0 = (2L«qJ[B] T [D][B]{q}-2L6qj[B] T [D]{e°} 

+ 2L6qJ[B] T {a°})dn-6LqJ{Q}. (A. 37) 

Rearranging (A. 37), 

[«q][/ fi ([B] T [D][B]{q}-[B] T [D]{e°}+[B] T {o 0 })dn-{Q}] = 0. (A. 38) 


122 


Since the [6q] are arbitrary variations, the quantity Inside the 
square brackets must vanish. Hence 

f Q ([B] T [D][B]{q})dn - {Q} + / n [B] T [D]{e°}dn 

- / n [B] T {o°}dn. (A. 39) 

Defining the element stiffness matrix [k] by 

[k] - [B] T [D][B]dn, (A. 40) 

(A. 39) becomes 

[k]{q> - {Q} + f a [B] T [D]{E°)dn - / Q [B] T {o°}dfi. (A.41) 

The last two terms In (A.41) represent the pseudo-loads necessary to 
account for the effects of Initial strains and Initial stresses. The 
equations (A.41) are the linear load-displacement relations for a single 
element. For an assemblage of elements, the equations (A.41) are 
"assembled" by superposition. This technique Is detailed In the text 
by Bathe and Wilson [63]. 

The integrations required in the evaluation of (A.41) are computed 
by Gaussian quadratures using a 4x4x2 rule. The Integrations to 
compute the element stiffness matrix will be discussed as an example. 
Gaussian quadrature Is discussed In detail in the text by 
Scarborough [77]. Considering again equation (A. 40), 

[k] - f Q [B] T [D][B]dft. (A. 40) 

When mapped Into the local €>n,c coordinate system, 
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dn - |J|dedndc (A.42) 

where |J| Is the determinant of the Jacobian [J]. The limits of the 
Integration are -1 to +1 In all directions. Equation (A. 40) becomes 

[k] - [B] T [D][B]|J|dcdndc (A.43) 

■ /It GU.n.Odedndc. 

The coordinates (a^.bj.c^) of the 32 Gauss points are 

a i - ±.86113, ±.33998 

bj * ±.86113, ±.33998 (A.44) 

c k - ±.57735 

The corresponding weighting functions ( H -j »Hj ,H k ) are 

H i - .34785 (for ±.86113), .65214 (tor ±.33998) 

Hj * .34785 (for ±.86113), .65214 (for ±.33998) (A. 45) 

H k « 1.0 (for ±.57735), 

Appling the rules of Gaussian quadratures, 

4 4 2 

[k] * l t Z H.H,H.G(a i ,b,,cJ (A.46) 

1=1 j«l k«l 1 J * 1 J K 

where G(e,n ,;) * [B] T [D][B] | J | . 

We note at this point that [k] Is a 72x72 matrix which Is symmetric. 

The other volume integral: required in (A. 41) are similarly evaluated. 
The result of the assembly process is a system of n linear 
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algebraic equations, where n Is the total number of degrees of 
freedom of the finite element assemblage. Before Imposition of 
boundary conditions sufficient to prevent rigid body motion, the 
system Is Indeterminate. After Imposition of boundary conditions, 
the global stiffness matrix [K] is synmetrlc, banded, and positive 
definite. Once all the terms In (A. 41) are evaluated, the nodal 
displacement vector {q } for the linear elastic problem is computed 
by Inversion of the global stiffness matrix. Thus, for linear elasto- 
s tat 1c problems, 

{q> - [k] _1 ({Q} + / n [B] T [D]{e°}dn - [B] T {o°)dn). (A.47) 

In practice, the global stiffness matrix is not actually Inverted. The 
system of equations is generally solved by some trlangularlzation and 
backsubstltutlon process. Bathe and Wilson [63] present a review of 
several commonly used equation solution techniques. 

For a nonlinear analysis using an Incremental intlal strain 
technique, Equ. (A.41), in the absence of Initial stresses, becomes 
for the !' t ^ 1 Increment 

[K]Uq} k - {AQ> k (A. 48) 

where 

{AQ) k = UQ} k ♦ {AQ 0 } k + { AQq } k + UQP} k ' f + (A. 49) 

and 

[K] ■ stiffness matrix 
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If 

UQ} * applied force Increment 

UQj) k ■ Incremental temperature pseudoload 

UQj) k ■ Incremental moisture pseudoload 

UQ^” 1 » Incremental plasticity pseudoload from prior step 

LI 

(Rq) « total equilibrium residual at end of prior step. 

The Incremental plasticity pseudoload Is designated (k-1) since It 
Is computed (see Chapter 3) using the plastic strain Increment computed 
during the (k-l) th Increment. The residual Is computed from a total 
equilibrium check at the end of t.ie (k-l) th Increment, and Is given 
by 

. , k-1 , k-1 , 

(Rr" 1 • [K]( E Uq} 1 ) - E UQ} 1 (A. 50) 

1-1 1-1 

and represents a correction which compensates for any errors In total 
equilibrium which may have existed at the end of step (k-1). The 
mechanics of computing the various pseudoloads are discussed In detail 
In Chapter 3. 
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APPENDIX B 

COMPUTER IMPLEMENTATION 

A computer program was written to test the theoretical developments 
previously presented. The program is called Nonlinear Analysis of 
Laminated Composites (NALCOM). This chapter presents a description of 
the cap > Titles of NALCOM as well as serves as a user's manual. 
Considerable effort was expended in choosing mnemonics consistent with 
the theoretical developments. The program is as modular as was deemed 
practical. Program development was done on the CDC 6700 computer at 
the United States Naval Surface Weapons Center, Dahlgren, Virginia. 

The program was later converted and run on the IBM S370/3303 at Virginia 
Polytechnic Institute and State University, Blacksburg, Virginia. Due 
to core limitations of the CDC system, considerable equivalencing is 
used, many parameters are passed in subroutine calls, and scratch files 
are used for storage of matrices, material properties, etc. 

B.l Summary of Capabilities 

The only finite element currently implemented in NALCOM is the 24 
node isoparametric element (Figure A.l) with transversely isotropic 
material properties. All capabilities are available for transversely 
isotropic materials. The nonlinear load-deflection equations are 
solved incrementally. Any combination of the three load types (with 
the exception of specified displacements with plasticity) may be 
applied to the finite element model, but only one type of load may 
be applied per load step. As with any nonlinear analysis program. 
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the choice of load step Is left to the analyst. Prospective users are 
encouraged to read and understand the theoretical basis of NALCOM. 
Violation of assumptions by an Improper choice of load step will result 
In erroneous results. NALCOM was developed as a research tool, and 
while some error checking Is done, other operations are not monitored 
for quality. To paraphrase, "If you don't know how It works, don't 
fool with It!" 

B.2 Program Description 

NALCOM Is written In ANSI standard FORTRAN and consists of 
approximately 4000 source statements, may of which are comments. A 
listing of the program Is available from the author. Core requirement 
for in-core solution on a CDC 6000 series computer Is approximately 
100000 octal words plus the number of stiffness matrix elements stored. 
The number of elements (or d.o.f.) which may be used In the finite 
element model Is based on the available storage, number of element 
groups (see input description), and the bandwidth of the global 
elastic stiffness matrix. The author has used a modified version of 
the NASTRAN preprocessor BANDIT [64] to minimize bandwith prior to 
any runs on the problems reported. Due to the use of the skyline 
column storage algorithm of Bathe and Wilson [63], minimization of 
mean bandwith is Important, since smaller mean bandwldths require less 
core storage and fewer arithmetic operations for trlangularlzatlon and 
back-substitution. 

NALCOM is specifically designed for nonlinear analysis. 
Consequently, linear analyses may be somewhat inefficient. For 


129 


example, cooling of a linear elastic material with constant properties 
followed by load application would ordinarily be handled by superposi- 
tion In one calculation. NALCOM, however, would require two load steps, 
first a cooling step, followed by a loading step. Such Inefficiency 
would not be acceptable In a production code, but Is believed acceptable 
In a research tool. Other Inefficiencies Include the carrying of all 
36 terms of elasticity and compliance matrices. Other quantities, such 
as the B matrices (strain-displacement transformation), determinant of 
the Jacobian for the Isoparametric mapping, stress and strain trans- 
formation matrices are stored in core or on scratch files to avoid 
recomputation. 

The starting point for the development of NALCOM was the sample 
program STAP given in the text by Bathe and Wilson [65]. The concept 
of element groups, with properties stored on disk when not in use; the 
skyline column storage scheme utilizing the efficient equation solver 
COLSOL; and the semi-dynamic storage allocation technique all made STAP 
a perfect starting point. For details on the skyline column storage 
scheme and the routine, the reader is directed to the text by Bathe 
and Wilson [63]. The program STAP was modified to use the 24 node 
isoparametric element. Additional modifications were made, in a step- 
by-step manner, to introduce all capabilities. Checkout has included 
a wide variety of problems, some of which are discussed in Chapter 5. 
Complete checkout of all options for all possible loading combinations, 
material models, and geometries is beyond the scope of the present 
study. It is possible that future NALCOM users will encounter some 
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difficulties. This possibility places additional emphasis on the need 
for understanding of the theoretical basis of the program. 

The high speed storage allocation for NALCOM is shown In Figures 
B.l. Note the similarity to Figure 6.6 of Reference 63. The 
significance of the arrays which are Input is noted In the Input guide 
presented later in this Appendix. For detailed discussions of the 
element connectivity array LM, the skyline diagonal pointer array MAXA, 
and the column height vector MHT, the reader is directed to the text 
by Bathe and Wilson [63]. The variable IT’JO is an aid to machine-to- 
machine conversion. It has a value of one for single precision and 
two for double precision. A flow chart for NALCOM is shown in 
Figure B.2. Table B.2 is a list, along with a brief description, 
of the scratch files used. 

Several comments are necessary regarding the implementation of 
the temperature dependent material properties and plasticity. All 
temperature dependency relationships are input as piecewise linear 
with a maximum of ten points (nine linear segments). The maximum 
of ten is easily changed by changing dimensions and loop ranges in 
the appropriate reading and interpolating subroutines. A minimum of 
two points is required, and if temperature dependent elastic constants 
are declared (ITL0AD.GE.2) , all properties are automatically tempera- 
ture dependent. Properties at points between the extremes are linearly 
interpolated on the appropriate segment. A temperature lower than the 
lowest input temperature will result in properties evaluated at the 
lowest input. Similarly, a temperature higher than the highest input 
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Address In 
Blank Common 

Storage 
Requl red 

Array 

Nl=l 

3*NUMNP 

ID array 

N2 

NUMNPMTWO 

X-coordinate array 

N3 

NUMNP*ITWO 

Y-coordinate array 

N4 

NUMNP*ITWO 

Z-coordinate array 

N5 




FIGURE B-la. INPUT PHASE STORAGE ALLOCATION 
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Address in Storage 

Blank Common Required 


Nl-1 

N2 

N3 


N4 


3*NUMNP 

NEQ+1 

NWK*ITW0 


NEQ*ITW0 


N5 MAXEST 

N6 


FIGURE B-lb. K ASSEMBLY, SOLUTION, STRESS 
ALLOCATION 


Array 


ID array 


MAXA array 


Global stiffness 
matrix K 


R-load vector or 
DU- displacement 
vector 


Element group 
Information 


COMPUTATION STORAGE 
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Address In 
Blank Common 

Storage 

Required 

Array 

Nl-1 

3*NUMNP 

ID array 

N2 

NEQ+1 

MAXA array 

N3 

NWK*ITWO 

Global stiffness 
matrix K 

N4 

NEQ*ITWO 

R-load vector 

N5 

MAXEST 

Element group 
information 

N6 

NLOAD 

NOD- loaded nodes 

N7 

NLOAD 

IDIRN- load 
directions 

N8 

NLOAD* I TWO 

FLOAD- load 
magnitudes 


N9 


FIGURE B-lc. LOAD VECTOR GENERATION PHASE STORAGE 
ALLOCATION 
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Address In 

Storage 

Array 

Blank Common 

Required 


Nl a l 

3*NUMNP 

ID array 

N2 

NEQ+1 

MAXA array 

N3 

NWK*ITWO 

Global stiffness 
matrix K 

N4 

NEQ*ITWO 

R-load vector 

N5 

MAXEST 

Element group 
Information 

N6 

NBC 

NOD- nodes with 
specified displace- 
ments 

N7 

NBC 

IDIRN-dlrectlon of 
specified displace- 
ments 

N8 

NBC*ITWO 

BC- boundary condi- 
tion vector 


N9 


FIGURE B-ld. SPECIFIED DISPLACEMENT DEFINITION PHASE 
STORAGE ALLOCATION 
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Address In Storage Array 

Blank Common Required 


Nl*l 

3*NUMNP 

ID array 

N2 

NEQ+1 

MAXA array 

N3 

NWK*ITWO 

Global stiffness 
matrix K 

N4 

NEQ*ITWO 

R-load vector 

N5 

NEQ*ITWO 

BC- boundary condi- 
tion vector 

N6 




FIGURE B-le. SPECIFIED DISPLACEMENT IMPOSITION PHASE 
STORAGE ALLOCATION 


136 


Address In 
Blank Common 

Storage 
Requl red 

Array 

N101 

5*ITW0 

E- room temperature 
elastic constants 

N102 

6*ITW0 

ALP- room temp thermal 
expansion coefficients 

N103 

6*ITW0 

not currently used 

N104 

I TWO 

ANG- ply orientation, 
degrees 

N105 

36* I TWO 

D- room temp elasticity, 
global coordinates 

N106 

36*ITW0 

012- room temp 
elasticity matrix, 
material coordinates 

N107 

36* I TWO 

COM- room temp compli- 
ance matrix, global 
coordinates 

N108 

36* I TWO 

C0M12- room temp compli- 
ance matrix, material 
coordinates 

NT 09 

36* I TWO 

Tl- stress transforma- 
tion matrix 

NT 1 0 

36* I TWO 

Til- Inverse of Tl 

Nlll 

36* I TWO 

T2- strain transforma- 
tion matrix 

Nil 2 

36* I TWO 

T2I- Inverse of T2 

Nil 3 

72*NUME*ITW0 

EGXYZ- element group 
coordinate array 

N114 

72*NUME 

LM- element connectivity 
array 

Nil 5 

NEQ*ITWO 

V- scratch area 

N116 




FIGURE B-lf. ELEMENT GROUP INFORMATION STORAGE ALLOCATION 
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READ CONTROL INFORMATION 





[ 

READ NOOAL COORDINATES, 
ESTABLISH EQUATION NUMBERS 




GENERATE GAUSSIAN QUADRATURE 
FACTORS, TRANSFORMATION MATRICES 





1 RE AO, STORE ELEMENT DATA 


I 


ASSEMBLE GLOBAL ELASTIC 
STIFFNESS MATRIX 


'read LOAD TYPE. NUMBER OF STEPS | 


| GENERATE. STORE LOAD VECTOR! 


1 

1 

) 

\ 

IMPOSE SPECIFIED DISPLACEMENTS 

IF REQUIRED 


I 


REREAD AND/OR TRIANGULARIZE 
STIFFNESS MATRIX 


[COMPUTE DISPLACEMENTS ) 
COMPUTE STRESSES, STRAINS 



FIGURE B.2. FLOW CHART FOR NALCOM 
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TABLE B-1 . NALCOM SCRATCH FILE ASSIGNMENTS 


File Name 

Unit Number 

Use 

IELMNT 

1 

element group Information 

I LOAD 

2 

Incremental load vector 

IBCTAP 

3 

specified displacement vector 

ICTAPE 

4 

coordinate transformation matrices 

I IN 

5 

Input device 

I OUT 

6 

output device 

ISTRS 

8 or 18 

stress Information 

ISTRN 

9 or 19 

strain Information 

ISTF 

14 

global elastic stiffness matrix 

IALPH 

15 

thermal expansion coefficients 

IPCR 

16 

percent retention data 

IRAMOS 

17 

Ramberg-Osgood coefficients 

ISTRSO 

8 or 18 

stress information 

ISTRNO 

9 or 19 

strain Information 

IDISP 

20 

total displacement vector 

IPLAS 

21 or 22 

plasticity data 

IPLASO 

21 or 22 

plasticity data 

IBTAPE 

23 

B (strain-displacement transformation) 



and DETJ (Jacobian determinant) data 

ISTFT 

24 

total mechanical force vector 

IPLTAP 

25 

plasticity pseudoload vector 
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will result In properties at the highest Input. 

Linear Interpolation Is also used to compute the fraction of the 
step which was elastic for an element undergoing the transition from 
elastic to plastic behavior. The Interpolation Is done on the numerical 
value of Hill's yield condition evaluated at the beginning and end of 
the load step. Since Hill's criterion Is quadratic In stress, this 
Interpolation Is accurate only for relatively small Increments. Also, 
too many Gauss points should not be allowed to transition during a 
single load step due to mutual Influence of the Gauss points on their 
neighbors. 

B.3 Program Input 

This section details the Inputs to NALCOM. It should be noted 
that no preprocessing Is currently Implemented. This Is due largely to 
core limitations of the COC 6700 on which the program was developed. 
Also, due to the previously mentioned Importance of bandwidth. It Is 
recommended that grids should be generated by some other means and the 
bandwidth minimized before use of NALCOM. 

B.3.1 Units 

There Is no dependence on the system of units used In NALCOM. Any 
set of consistent units may be used, but care must be taken that 
consistency Is present. The one exception to this rule Is that the 
ply orientation angle, ANG, must be Input In degrees. 

B.3. 2 Inputs, Formats, Mnemonics, and Description 
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I. TITLE CARO (20A4) 

Note Columns Variable Description 

(1) 1-80 n£0(20) 80 columns of alphanumeric title Info 


NOTES - 

(1) Two blank cards must be Input at the end of the data deck for 
normal termination 
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II. MASTER CONTROL CARD (715) 


Note 

Col umns 

Variable 

Description 

(1) 

1-5 

NUMNP 

Total number of nodal points, .EQ.O for 
program stop 

(2) 

6-10 

NUMEG 

Number of element grotps, .GT.O 

(3) 

11-15 

NLCASE 

Number of load cases, .GT.O 


16-20 

M0DEX 

Solution mode ; 

. EQ. 1 execution, no plasticity 
.EQ.2 execution, with plasticity 


21-25 

IPRINT ( 1 ) 

Equation number print flag, 

.EQ.O print equation numbers corres- 
ponding to nodal d.o.f. 

.EQ.l no print equation numbers 


26-30 

IPRINT ( 2 ) 

Gauss point print flag 

.EQ.O Compute and print Cartesian coord- 
inates of Gauss points 
.EQ.l No compute 

(4) 

31-35 

NITER 

Maximum number of equilibrium iterations 
per increment. Set .EQ.O for no iteration 


NOTES - 


(1) Controls number of cards of Type IV to be read 

(2) An element group is a convenient collection of elements, all 
of which must have the same fiber orientation for orthotropic 
materials 

(3) Number of unique load types or magnitudes to be applied 

(4) Convergence occurs when the change in the Euclidea.' norm 
of the residual from the last iteration is less than 0.001 
times the Euclidean norm of the incremental load vector 


I 
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III. STRESS FREE TEMPERATURE (F10.0) 
Note Columns Variable 


t 

t 




Description 


1-10 


SFT 


Laminate stress free temperature or ref- 
erence temperature for thermal expansion 




Note 


"1 
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IV. INPUT SCALE FACTORS (3F10.0) 


Col umns 

Variable 

Description 

1-10 

XFAC 

Input x coordinates will be multiplied 
by this quantity 

11-20 

YFAC 

Input y coordinates will be multiplied 
by this quantity 


21-30 


ZFAC 


Input z coordinates will be multiplied 
by this quantity 
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V. NODAL POINT CARDS (4I5.3F10.0) 


Note 

Columns 

Variable 


Description 

(1) 

1-5 

N 

Nodal point number, .GE.O and .LE.NUMNP 

(2) 

6-10 

ID(1,N) 

X 

translation code 


11-15 

ID(2,N) 

Y 

translation code 


16-20 

ID(3,N) 

Z 

translation code 


21-30 

X(N) 

X 

coordinate 


31-40 

Y(N) 

Y 

coordinate 


41-50 

Z(N) 

Z 

coordinate 


NOTES - 

(1) Data must be defined for all nodal points. Nodes may be 
input in any order. 

(2) Translation codes define whether a node is free to move. 
ID(M,N)=0 means node N may translate in the M direction 
ID(M,N)=1 means node N is fixed in the M direction. 

Equation numbers are not assigned to degrees of freedom with 
ID=1 . DOF with ID=1 are used to define rigid boundaries, 
symmetry lines, etc. Note that sufficient fixities must be 
present to prevent rigid body motion of the structure. 

(3) The XYZ coordinate system is hereafter referred to as the 
global, or laminate coordinate system. 
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VI. ELEMENT GROUP DATA 

The following sequence of cards Is present for each of the NUMEG 
element groups. The only element currently implemented is the 24 node 
Isoparametric element. Note that each element group may be a different 
material, different thermal, plasticity behavior, etc. 

VI.l Element Group Control Card (515) 


Note 

Col umns 

Variable 


1-5 

NPAR(l) 

(1) 

6-10 

NPAR(2) 
or NUME 


11-15 

NPAR(3) 
or ITLOAD 


16-20 

21-25 NPAR(5) 

or I DEN 


Description 


Enter the number 1 

Number of elements in this group .GE.l 


Thermal expansion type 
.EQ.O No thermal effects 
.EQ.l Constant coefficients of thermal 
expansion 

. EQ. 2 First order nonlinear thermal 
effects 

. EQ. 3 Second order nonlinear plus first 
order nonlinear thermal effects 

Not currently in use 

Element group identity flag, 

.EQ.O elements may be geometrically 
different 

.EQ.l elements are geometrically identical 


NOTES - 

(1) Element numbers begin with 1 and end with NUME within each 
element group. Tnere is no limit, other than available core, 
on the value of NUME. 
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VI. 2 ELASTIC PROPERTIES AND PLY ORIENTATION (6F10.0) 


Note 

Col umns 

Variable 

Description 

(1) 

1-10 

E(l) 

Longitudinal Young's modulus, E^ 


11-20 

E(2) 

Transverse Young's modulus, E 22 


21-30 

E(3) 

Shear modulus, 


31-40 

E(4) 

Major Poisson's ratio, v^ 2 


41-50 

E(5) 

Minor Poisson's ratio, 


51-60 

ANG 

Ply orientation (angle) measured from 
X- (global) to 1- (longitudinal material) 
axis, DEGREES (see Figure B.3) 


NOTES - 


(1) Properties are input as room temperature (RT) values 
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VI. 3 RAMBERG-OSGOOD COEFFICIENTS 

Two forms of Ramberg-Osgood coefficient Inputs are available, the 
first constant with respect to temperature (ITL0AD.LE.1) and the other 
variable with respect to temperature (ITL0AD.GE.2). Recall the 
Ramberg-Osgood form Is e = | + 8o n . 

Ramberg-Osgood coefficients are read only If M0DEX.EQ.2. 

VI. 3.1 Constant Ramberg-Osgood Coefficients (6E13.5) 

These cards are read If ITL0AD.EQ.1. Three cards are required, 
the first for the material longitudinal direction, the second for the j 

transverse direction, and the third for shear. Transition from n-j j 

and B.j to n 2 and e % occurs at o*. I 


Note 

Column 

Variable 

Description 

(1) 

1-13 

ROB(I) 

Ramberg-Osgood coefficient 8-j 


14-26 

RON(I) 

Ramberg-Osgood coefficient n^ 


27-39 

YIE(I) 

Yield stress In I direction 


40-52 

R0B( 1+4) 

Ramberg-Osgood coefficient Bg 


53-65 

RON (1+4) 

Ramberg-Osgood coefficient n 2 


66-78 

YIE( 1+4) 

Transition stress, a* 


NOTES - 


(1) Directions are given by 


1*1 longitudinal direction 
1=2 transverse direction 
1*3 shear direction 
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VI. 3. 2 Temperature Dependent Ramberg-Osgood Coefficients 
These cards are read If ITL0AD.GE.1. Temperature dependent 
Ramberg-Osgood coefficients and yield strengths are Input at up to ten 
temperatures, as specified on the following card. 

VI. 3. 2.1 Number of Temperatures for Ramberg-Osgood Input (15) 


Note 

Col umn 

Variable 

Description 


1-5 

NTEMS 

Number of temperatures for R-0 and yield 
strength Input, .GE.2 and .LE.10 

NTEMS sets of the four R-0 and yield data cards are required. In order 

of ascending temperature. 



VI. 3. 2. 2 

Temperature for R-0 and Yield Strength Input (F10.0) 

Note 

Col umn 

Variable 

Description 

(1) 

1-10 

TEMP 

Temperature 


VI. 3. 2. 3 

R-0 and Yield Data (6E13.5) 

Note 

Column 

Variable 

Description 

(1) 

1-13 

ROB(I) 

$1 at TEMP 


14-26 

RON(I) 

n ] at TEMP 


27-39 

YIE(I) 

Yield strength at TEMP 


40-52 

R0B( 1+4) 

B 2 at TEMP 


53-65 

RON (1+4) 

n 2 at TEMP 


66-78 

YIE(I+4) 

o* at TEMP 


NOTES - 


(1) Three of these cards are required at each value of TEMP, one 
for longitudinal, transverse, and shear behavior, In that 
order. If a bilinear R-0 fit is not desired, leave Columns 
40-78 blank. 
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VI. 4 THERMAL EXPANSION COEFFICIENTS 

Thermal expansion coefficients are read only If ITLOAD.GT.O. If 
ITL0AD.EQ.1, properties are constant. If ITL0AD.GE.2, properties vary 
with temperature, and piecewise linear (up to ten points) curves are 
Input. 

VI. 4.1 Constant Thermal Expansion Coefficients (2F10.0) 

Read If ITL0AD.EQ.1 . 


Note 

Column 

Variable 

Description 


1-10 

ALP(l) 

Longitudinal thermal expansion coefficient 


11-20 

ALP ( 2 ) 
=ALP(3) 

Transverse thermal expansion coefficient 


VI. 4. 2 Temperature Dependent Thermal Expansion Coefficients 

Read 

If ITL0AD.GT.1 . 



VI. 4. 2.1 

Number of Temperatures for Thermal Expansion Coefficient 
Input (15) 

Note 

Col umn 

Variable 

Description 


1-5 

NTEMS 

Number o^ temperatures for thermal expan- 
sion Input, .GE.2 and .LE.10 


VI. 4. 2. 2 

Variable Thermal Expansion Coefficient Cards (3F10.0) 
NTEMS cards are required. In order of ascending 
temperature. 

Note 

Column 

Variable 

Description 


1-10 

ALP ( 1 ) 

Longitudinal thermal expansion coefficient 
at TEMP 


11-20 

ALP(2) 

S ALP(3) 

Transverse thermal expansion coefficient 
at TEMP 


21-30 

TEMP 

Temperature 
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VI. 5 FRACTION RETENTION OF ROOM TEMPERATURE PROPERTIES 
If ITLOAD.GE.2, elastic coefficients are assumed to be variable 
with temperature. This Information Is Input as the fraction of room 
temperature (RT) properties retained at a given temperature. Piecewise 
linear curves of up to ten points are used. 

VI. 5.1 Number of Temperatures for Fraction Retention Input (15) 


Note Col umn 

Variable 

Description 

1-5 

NTEMS 

Number of temperatures for fraction 
retention Input, .GE.2 and .LE.10 

VI. 5. 2 Fraction Retention Data (6F10.0) 

Note Column 

Variable 

Description 

1-10 

TEMP 

Temperature 

11-20 

PEI 

Fraction retention of E^ at TEMP 

21-30 

PE2 

Fraction retention of E 2 2 at TEMP 

31-40 

PG12 

Fraction retention of G^ at TEMP 

41-50 

PGNU12 

Fraction retention of at TEMP 

51-60 

PGNU23 

Fraction retention of v 23 at TEMP 
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VI. 6 ELEMENT CONNECTIVITY INFORMATION 

Two cards are required for each element. Elements are numbered 
from ONE (1) to NUME within each element group* and must be Input In 
Increasing order. 


VI. 6.1 Element Card 1 (1315) 

Note Column Variable Description 


1-5 

6-10 

n-i5 

16-20 

etc. 

N 

NCON(I) 
1*1 ,12 

Element number, .GE.l and .LE.NUME 
First 12 node numbers of element N 

(See Figure A.l for Input order) 

VI. 6. 2 Element Card 2 (5X.12I5) 

Note Col umn 

Variable 

Description 

1-5 

— 

Leave blank 

6-10 

NCON(I) 

1*13,24 

Last 12 node numbers for element N 
(See Figure A.l) 
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VII. LOAD CASE DATA 

Data from this section must be present for each of the NLCASE load 
cases. Note that only one type of load Increment (force, temperature 
change, or specified displacement) may be applied In a load case. Also 
note thut specified displacements are Incompatible with plasticity as 
Implemented In NALCOM. 


VI 1.1 Load Case Control 

(315) 

Note Column 

Variable 

Description 

1-5 

LDTYP 

Load type, 

■1 , concentrated force 
>2, temperature change 
>4, specified displacement 

6-10 

NSTP 

Number of equal steps of this load type 
to be applied in this load case 

11-15 

IPRINT ( 3 ) 

Print flag, 

*0, print displacement, strain, and 
stress Information for every load 
step of this load case 
■1 , only print Information for last load 
step of this load case 
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VII. 2 Concentrated Force Load (LDTYP.EQ.l) 

VII. 2.1 Number of Loaded Degrees of Freedom (15) 


Note 

Column 

Variable 

Description 


1-5 

NLOAD 

Number of loaded degrees of freedom 


VII. 2. 2 

Concentrated Force Card(s) (2I5.F10.0) 
NLOAD of these cards are required. 

Note 

Column 

Variable 

Description 


1-5 

NOD 

Node to which force Is applied 


6-10 

IDIRN 

Direction in which force acts, 
■1, X-dlrectlon 
e 2, Y-directlon 
*3, Z-dlrectlon 


11-20 

FLOAD 

Increment of load to be applied for 
each of the NSTP steps of this load 
case 


VII. 3 Temperature Increment (LDTYP.EQ.2) (F10.0) 

Note 

Column 

Variable 

Description 

(1) 

1-10 

DTEMP 

Temperature Increment to be applied for 
each of the NSTP steps of this load case 


NOTE - 

(1) Only uniform temperatures are currently Implemented, so all 
temperatures are constant through the entire structure. 
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VII. 4 Specified Displacement Loads (LDTYP.EQ.4) 

VI. 4.1 Number of Degrees of Freedom with Specified 
Displacements (15) 


Note 

Col umn 

Variable 

Description 


1-5 

NBC 

Number of degrees of freedom with 
specified displacements 


VII. 4.2 

Specified Displacement Card (215 .FI 0.0) 
NBC of these cards are required. 

Note 

Col umn 

Variable 

Description 

(1) 

1-5 

NOD 

Node to which specified displacement is 
applied 

(2) 

6-10 

IDIRN 

Direction of displacement 
=1, X-direction 
=2, Y-dlrection 
=3, Z-direction 

(3) 

11-20 

DISP 

Incremental displacement to be applied for 
each of the NSTP steps of this load case 

NOTES 

. 




(1) Specified displacement loads are not compatible with 
plasticity 

(2) ID( IDIRN »N0D) must be nonzero, i.e. an equation must exist 
for this degree of freedom 

(3) The value of DISP may be zero 
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B.4 Program Output 

Printed output consists of an echo of the input data as it Is 
read, Cartesian coordinates of the Gauss points (for IPRINT(2).EQ.O), 
stiffness matrix bandwidth information, nodal displacements, and 
stress and strain data at the Gauss points. The amount of output is 
controlled by the values in the IPRINT array. Strain output consists 
of total strains at the Gauss points of each element. Strains are 
output in the global (XYZ) and material (123) coordinate systems. 
Plasticity output includes the global plastic strain increment, value 
of Hill's yield function (YF) at the end of the step, the maximum 
value Hill's yield condition has achieved (YFMAX), the fraction of 
the load step which was elastic (PCEL), and the loading parameter 
(LAMBDA). Also printed is an indicator (NYC), which has the following 
meanings: 

NYC = -1 Gauss point unloading elastically or reloading 
elastically toward subsequent yield 
NYC = 0 Gauss point is, and has always been, elastic 

NYC = 1 New (first step) yield or subsequent yield. Part 

of the stpp may have been elastic 
NYC = 2 Prior plastic, continuing plastic flow. 

The dominant direction of the plastic flow is denoted by the variable 
I DOM. 

It should be noted that no postprocessing is currently available. 
Scratch files IDISP, ISTRS, and ISTRN, as well as others, are available 
at the end of the run, and could be read by a postprocessor. A 
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technique of this type was used for generation of some of the plots 
In the body of the report. 


r 
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APPENDIX C 

YIELD AND PLASTICITY RELATIONS FOR ISOTROPIC MATERIALS 
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APPENDIX C 

YIELD AND PLASTICITY RELATIONS FOR ISOTROPIC MATERIALS 


The purpose of this appendix Is to show how Hill's orthotropic 
yield function (equs. 3.22 and 3.23), equivalent stress (equ. 3.24), 
and equivalent plastic strain increment (equ. 3.27) reduce to the 
familiar von Mises relationships for Isotropic materials. If these 
can be shown, the relationships 3.56-3.60 follow, since the same 
heuristic arguments made for orthotropic materials can also be 
advanced for isotropic materials. 

Consider Hill's orthotropic yield function, 

f( Si ) = F(S 2 -S 3 ) 2 + G(S 3 -S 1 ) 2 + H(S r S 2 ) 2 

+ 2LS 2 3 + 2MS 2 3 + 2NS 2 2 =1, (C.l 


or _ 1 . 1 1 

ZF " T + T " T 

Y* V- r 

?r _ 1 . 1 1 

V V 

2H = — i + — i- - ~ y 

r i c 


2l = T 
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and X, Y, Z, R, S, and T are the yield strengths in the principal 


material directions 1 

, 2, 3, 23, 13, and 12, respectively. 

For an 

isotropic material, 

X = Y = Z 

(C.3) 

and 

R - S = T. 

(C.4) 

Also, Mendel son [38] 

shows that for a distortion energy theory of 

yielding. 

R = — . 

(C.5) 


X 

Therefore, 

F = G « H » -!y 
2k 

(C.6) 

and 

2L = 2M = 2N = 4 

(C.7) 


Substituting C.6 and C.7 into C.l results in 


f < s i> ’ ^< S 2- S 3> 2 * < S 3- S 1> 2 + (S r S 2 )2] 

* ^ S Z3 tS 13 +S 12^ ’ 

A 

thus 

^(s 1 ) - C(S 2 -S 3 )2 + (Sj-S,) 2 + (S } -S 2 ) Z 

X 

+ 6s| 3 + 6S* 3 + 6S^ 2 ] = 1, (C.9) 
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where X Is the yield stress In simple tension. Note that C.9 is the 
familiar von Hises yield criterion for an isotropic material. 

Next consider the effective stress, 5 2 of equ. 3.24, 

52 - r f < s 2- s 3> 2 ♦ g < s 3-V 2 + H < s rV 2 


+ 2LS| 3 + 2MS^3 + 2NS* 2 ]. 

Note that 5 2 Is related to f(s^) by 

52 a 2{M+h') f(s 1 ) ‘ 

Using C.6 and C.9 In conjunction with C.ll yields 

52 = I C ] ^ t( S 2 _S 3 )2 + (S 3' S 1 )2 + ( S 1 _S 2 )2 


(C.10) 


(C.ll) 


and 


+ 6S 


23 


6S 13 + 6S 12^’ 


(C.12) 


5^ = 


(S 2 -S 3 ) 2 + (S 3 -S 1 ) 2 + (S r S 2 ) 2 + 6S 2 3 + 6S 2 3 + 6S 2 2 - 


(C.13) 


Equation C.13 is the von Mlses or effective stress for an isotropic 
material. Finally, consider the equivalent plastic strain Increment 
of equ. 3.27, 


d§ P » | (F+G+H) 

, 6 < H<i «33- Fde n t 
(FG+GH+HF) 2 


F(Gde 22 -Hde5 3 ) 2 

(FG+GH+HF) 2 

H(Fde P r Gde P 2 ) 2 

(FG+GH+HF) 2 
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< de 23> 2 

5l 


(de'j) 2 ( de i2> 2 

+ -r B? — t- 5H — 


1 

7 


(C. 14) 


p 

where de^j tyj have been converted to engineering shear strain. 
Substituting C.6 and C.7 Into C.14 results In 


di p 


2 / 3 » 
6 IT 


% ^ 2 [< de 2 P r de «) 2 + KX/ * (de n' de 22 ): 

(-V 


+ T" C( de 23> 2 + ( de 13> 2 + ( d ^ 2 ) 2 ] 1 (c.15) 

d 5 P * l 2 -T XXs^ + ( d 4- de l P l> 2 + (de ir de 22 )2] 

* ^ C(de P 3 ) Z ♦ (de P 3 ) Z + (de P 2 ) 2 ] 7 (C.16) 


di P = [( de 22" de 33^ 2 + ( de 33 -de ll) 2 + ( de iT de 22^ 2 

1 

+ 6(de P 3 ) 2 + 6(de P 3 ) 2 + 6(de P 2 ) 2 ] ? (C.17) 

which Is the same form of the equivalent plastic strain Increment for 
an isotropic material, where the shear strains are engineering shear 
strains. 

The relationships between Hill's orthotropic plasticity model 
and the von Mises (or distortion energy) theory for Isotropic 
materials have been shown. The model developed In Chapter 3 can 
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be used for Isotropic as well as orthotropic materials. Results from 
Isotropic material analyses exhibit good experimental /analytical 
correlation, as shown In Figure 5.1. 
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